Hybrid Quantum Simulation

Hybrid quantum simulation deliberately composes quantum and classical representations, or distinct modes of quantum control, so that each component handles a part of a specified simulation task. The word hybrid describes the architecture of the workflow. It does not imply low depth, near-term utility, variational optimization, or freedom from error.

Three families account for most uses of the term:

1. Digital–analog simulation interleaves continuously parameterized native many-body evolution with discrete gates or pulses.
2. Variational quantum–classical simulation estimates a projected equation of motion or objective on a quantum device and updates a compact state representation classically.
3. Quantum embedding assigns a correlated fragment or impurity to a quantum solver while a classical environment supplies boundary data and enforces a closure or self-consistency condition.

These families can overlap. A variational impurity solver may itself use digital–analog control. A digital–analog state-preparation block may sit inside a larger classical feedback loop. Classification is useful only when it makes the interfaces, resources, and failure modes more explicit.

Required background. What Is Quantum Simulation? defines the target question; digital and analog simulation explain its two main implementations. Unitary evolution supplies the propagators used below.

Each approach must preserve the meaning of exchanged quantities: a pulse angle needs its rotation convention, a Green’s function its orbital ordering and frequency units, and several estimated observables their covariance. A shallow circuit can save coherent evolution time while requiring many adaptive measurements. A converged embedding calculation can still contain environment-model error. The relevant comparison is the accuracy and total cost of the requested physical output.

Digital–Analog Quantum Simulation

A digital–analog simulator composes native many-body evolutions with discrete control operations. Let $H_j(\boldsymbol\lambda_j)$ be a calibrated native Hamiltonian available for duration $\tau_j$, and let $D_j$ denote a discrete gate layer or fast control pulse. A general sequence has the form

$$U_{\mathrm{DA}} = D_m e^{-iH_m\tau_m/\hbar} \cdots D_2 e^{-iH_2\tau_2/\hbar} D_1 e^{-iH_1\tau_1/\hbar}.$$

The continuously tunable entangling evolutions are the analog blocks; the discrete operations select frames, refocus terms, alter signs, permute degrees of freedom, or compose blocks into a broader model family. The sequence may be open loop. No classical optimizer is required for this use of hybrid.

Toggling frames

Conjugating an available block by a digital operation changes its generator:

$$D e^{-iH\tau/\hbar}D^\dagger = e^{-i(DHD^\dagger)\tau/\hbar}.$$

This identity is the basic mechanism behind sign reversal, basis changes, and Hamiltonian refocusing. If a hardware interaction supplies $ZZ$ couplings, single-qubit rotations can turn them into $XX$ or $YY$ couplings. Local $X$ pulses can reverse selected $Z$ fields or selected Ising bonds. The useful question is not whether the native block is universal by itself, but which target generators lie in the control algebra generated by the available blocks and discrete operations at acceptable cost.

Product composition and noncommutativity

For two bounded generators $A$ and $B$, a short block pair obeys the Baker–Campbell–Hausdorff expansion

$$\begin{aligned} e^{-iA\Delta t/\hbar} e^{-iB\Delta t/\hbar} &= \exp\!\left[ -\frac{i\Delta t}{\hbar}(A+B) -\frac{\Delta t^2}{2\hbar^2}[A,B] +O(\Delta t^3) \right]. \end{aligned}$$

Thus the effective generator is

$$H_{\mathrm{eff}} = A+B -\frac{i\Delta t}{2\hbar}[A,B] +O(\Delta t^2).$$

Here $H_{\mathrm{eff}}$ generates one target-time step $\Delta t$. If each native block runs for laboratory duration $\Delta t$, the pair occupies $2\Delta t$ before pulse and switching overhead. Its generator per unit laboratory time is therefore $H_{\mathrm{eff}}/2$. In the limit of many short pairs, simulating $A+B$ for time $t$ still takes laboratory time $2t$ under these assumptions; refining the steps does not remove this factor. Decoherence and heating accumulate over the actual execution time. The amplitude-damping benchmark shows how to convert a laboratory decay rate to the target clock.

Noncommuting blocks therefore introduce a controllable composition error much like a product formula. That algorithmic error is distinct from mismatch in the physical generators $A$ and $B$. Increasing the number of blocks can reduce the former while increasing exposure to pulse error, drift, leakage, and decoherence.

The abstract order conditions are derived in Trotter Product Formula, while Trotter–Suzuki Methods develops simulation-level commutator bounds, ordering, cost, and refinement tests. Here the important point is architectural: a nominally exact analog block does not make the composed sequence exact.

Worked example: an exact refocusing block

Suppose the native two-qubit Hamiltonian is

$$H = J Z_1Z_2 +h_1Z_1 +h_2Z_2.$$

Let $P=X_1X_2$. Because $XZX=-Z$,

$$PHP = JZ_1Z_2 -h_1Z_1 -h_2Z_2.$$

All terms are diagonal in the computational basis, so $H$ and $PHP$ commute. The echo sequence

$$U_{\mathrm{echo}}(\tau) = P e^{-iH\tau/(2\hbar)} P e^{-iH\tau/(2\hbar)}$$

is therefore exactly

$$\begin{aligned} U_{\mathrm{echo}}(\tau) &= e^{-i(PHP)\tau/(2\hbar)} e^{-iH\tau/(2\hbar)} \\ &= e^{-iJZ_1Z_2\tau/\hbar}. \end{aligned}$$

Under the ideal model, the local fields cancel with no product-formula error while the interaction survives. In a laboratory, finite pulse duration can allow $H$ to act during $P$; pulse-axis error can leave residual fields; extra $XX$, $ZX$, or spectator couplings may not commute; and the two half-intervals may see different drift. The exact algebra establishes the intended block, not its physical fidelity.

Stepwise and always-on control

In a stepwise protocol, analog interactions are nominally switched on and off between digital layers. In an always-on or banged protocol, fast digital pulses are applied while a background interaction continues to act. The latter can avoid switching transients, but the digital and analog controls then overlap. An instantaneous-pulse model must be replaced by a time-ordered Hamiltonian,

$$U(T) = \mathcal T \exp\!\left[ -\frac{i}{\hbar} \int_0^T \bigl(H_{\mathrm{native}}+H_{\mathrm{ctrl}}(t)\bigr) \,dt \right].$$

Whether one version is better is hardware- and task-dependent. Switching errors, bandwidth, pulse distortions, crosstalk, and analog-block calibration must enter the comparison alongside circuit depth.

Errors in composed native evolutions

An ideal block sequence can differ from the target through its synthesis approximation. Physical execution adds generator mismatch, finite-pulse effects, decoherence, leakage, and preparation/readout errors. These terms can be combined only after expressing them in compatible norms or final-observable units.

For ideal and implemented unitary blocks with operator-norm errors $\epsilon_j$, a telescoping argument gives the conservative bound

$$\left\| \widetilde U_m\cdots\widetilde U_1 -U_m\cdots U_1 \right\| \le \sum_{j=1}^{m}\epsilon_j.$$

This is not a prediction that every error adds linearly. Coherent errors can cancel or align, stochastic errors can accumulate differently, and observable- specific error can be much smaller. The bound does show why counting native blocks without characterizing them is not a resource or accuracy estimate.

Variational Quantum Simulation

Variational quantum simulation compresses a state or process into a parameterized manifold and uses measured quantities to choose its motion. It is different from merely optimizing a final energy: the classical update is intended to approximate a differential equation for real time, imaginary time, or an open-system process.

Let

$$\lvert\psi(\boldsymbol\theta)\rangle \in \mathcal M$$

be a normalized ansatz with real parameters $\theta_i$. Exact Schrödinger evolution generally points outside the tangent space of $\mathcal M$. A variational principle projects the exact velocity onto the accessible tangent directions, after removing the physically irrelevant global phase.

Phase-invariant tangent projection

Define the projective tangent vectors

$$\lvert D_i\psi\rangle = \left( I-\lvert\psi\rangle\langle\psi\rvert \right) \lvert\partial_i\psi\rangle,$$

and the instantaneous energy

$$E = \langle\psi\rvert H\lvert\psi\rangle.$$

McLachlan projection minimizes the norm of the projective residual

$$\lvert R\rangle = \sum_j \dot\theta_j \lvert D_j\psi\rangle +\frac{i}{\hbar} (H-E) \lvert\psi\rangle.$$

Stationarity with respect to each real velocity $\dot\theta_i$ gives

$$\sum_j A_{ij}\dot\theta_j = C_i,$$

with

$$\begin{aligned} A_{ij} &= \operatorname{Re} \langle D_i\psi\vert D_j\psi\rangle, \\ C_i &= \frac{1}{\hbar} \operatorname{Im} \langle D_i\psi\rvert H\lvert\psi\rangle. \end{aligned}$$

$A$ is the real part of the quantum geometric tensor on the ansatz manifold. It is positive semidefinite and can be singular when parameters are redundant or locally inactive. Some presentations fix a phase convention instead of using $\lvert D_i\psi\rangle$ explicitly; sign and gauge conventions must be kept consistent when formulas are combined.

At one time step, the hybrid loop is:

1. prepare $\lvert\psi(\boldsymbol\theta_k)\rangle$;
2. estimate the required entries of $A$ and $\boldsymbol C$;
3. solve a regularized linear system for $\dot{\boldsymbol\theta}_k$;
4. integrate to $\boldsymbol\theta_{k+1}$;
5. repeat, while evaluating validation observables and stopping conditions.

The general design of ansatzes, gradients, shot allocation, and optimizers is covered by Variational Quantum Algorithms. The distinctive issue here is whether a measured tangent projection tracks the target dynamics over time.

Projection error is a physical diagnostic

The minimum local residual

$$\epsilon_{\mathrm{proj}}(t) = \min_{\boldsymbol v} \left\| \sum_j v_j\lvert D_j\psi\rangle +\frac{i}{\hbar}(H-E)\lvert\psi\rangle \right\|$$

measures how much of the exact projective velocity lies outside the ansatz tangent space. A small optimization residual for the linear solve is not the same quantity. One can solve $A\dot{\boldsymbol\theta}=\boldsymbol C$ perfectly even when the ansatz manifold misses most of the target velocity.

The local projection residual is also not the entire trajectory error. Time integration, finite sampling, regularization, hardware bias, and repeated feedback can accumulate. Observable-specific checks along the path remain necessary.

Worked example: exact one-qubit motion

Take

$$H = \frac{\hbar\Omega}{2}Z, \qquad \lvert\psi(\theta)\rangle = e^{-i\theta Z/2} \lvert +x\rangle.$$

Because $\langle Z\rangle=0$ on this orbit,

$$\lvert D_\theta\psi\rangle = -\frac{i}{2}Z\lvert\psi\rangle.$$

The scalar metric and force are

$$A = \frac14, \qquad C = \frac{1}{\hbar} \operatorname{Im} \langle D_\theta\psi\rvert H\lvert\psi\rangle = \frac{\Omega}{4}.$$

Therefore $\dot\theta=\Omega$, exactly matching Schrödinger evolution. The ansatz is invariant under the target Hamiltonian, so $\epsilon_{\mathrm{proj}}=0$. For this constant-$\Omega$ parameter equation, the Euler update $\theta_{k+1}=\theta_k+\Omega\Delta t$ is itself exact; evaluating the ansatz at the updated angle gives every observable exactly in ideal arithmetic. Nonlinear parameterizations or state-dependent velocities can introduce integration error, while measured estimates of $A$ and $C$ still carry uncertainty.

Ill-conditioning amplifies estimator error

Suppose measurements return

$$\widehat A=A+\delta A, \qquad \widehat{\boldsymbol C} = \boldsymbol C+\delta\boldsymbol C.$$

On a fixed retained parameter subspace where $A$ is invertible, small perturbations give, to first order,

$$\delta\dot{\boldsymbol\theta} \approx A^+ \left( \delta\boldsymbol C -\delta A\,\dot{\boldsymbol\theta} \right),$$

where $A^+$ denotes the inverse on that subspace and zero on discarded directions. This formula does not include changes in a data-dependent retained subspace or rank; near a truncation threshold, those changes must be analyzed separately. Consequently,

$$\left\| \delta\dot{\boldsymbol\theta} \right\| \lesssim \left\|A^+\right\| \left( \left\|\delta\boldsymbol C\right\| + \left\|\delta A\right\| \left\|\dot{\boldsymbol\theta}\right\| \right).$$

Near-null tangent directions make $\lVert A^+\rVert$ large. Truncating small eigenvalues, adding Tikhonov regularization, or constraining the step can reduce variance, but each changes the projected dynamics and introduces bias. Metric eigenvalues also depend on parameter units and scaling. Check the induced state velocity and physical observables, not only large coordinate velocities. The threshold and coordinates are part of the numerical approximation.

Real time, imaginary time, and general processes

For normalized imaginary-time evolution with inverse-energy coordinate $\tau$,

$$\frac{d}{d\tau} \lvert\psi\rangle = -(H-E)\lvert\psi\rangle,$$

the same tangent projection yields

$$\sum_j A_{ij}\dot\theta_j = -\operatorname{Re} \langle D_i\psi\rvert H\lvert\psi\rangle.$$

Imaginary-time projection can prepare low-energy states, but finite ansatz capacity, normalization, conditioning, and measurement cost still control its quality. VQE is the derivation of static variational energy bounds.

Variational principles can also be formulated for density operators, purifications, stochastic trajectories, and general linear processes. In those settings, preserving trace, Hermiticity, positivity, and complete positivity may require more than minimizing a Hilbert-space residual. A low-residual parameter trajectory is not by itself a certificate that every intermediate object is a physical state or channel.

How local errors accumulate

Ansatz restriction, tangent projection, finite-shot estimates, regularization, time integration, and hardware error affect different parts of the trajectory. They are not independent additive errors: changing the ansatz changes the tangent space, while hardware bias changes the estimated vector field followed by the integrator. Step-size convergence should therefore be tested separately from ansatz enlargement and measurement precision.

Quantum Embedding and Self-Consistency

Quantum embedding assigns a selected correlated region to a quantum solver and represents the rest of the system through a classical environment. The full Hamiltonian may be partitioned schematically as

$$H = H_F +H_E +H_{FE},$$

where $F$ denotes a fragment or impurity and $E$ its environment. The method does not usually discard $E$. It compresses the environment into boundary data such as a bath hybridization, self-energy, correlation potential, mean field, or reduced density matrix.

An abstract embedding loop is

$$\begin{aligned} H_{\mathrm{imp}}^{(k)} &= \mathcal B\!\left[x_k\right], \\ \widehat y_k &= \mathcal Q\!\left[H_{\mathrm{imp}}^{(k)}\right], \\ x_{k+1} &= \mathcal F\!\left(x_k,\widehat y_k\right). \end{aligned}$$

$\mathcal B$ constructs the finite quantum problem from classical environment data $x_k$. The quantum solver $\mathcal Q$ returns observables $\widehat y_k$. The classical closure $\mathcal F$ updates the environment. At convergence, the selected fragment and environment quantities satisfy a prescribed matching condition.

DMFT and DMET illustrate different interfaces

In dynamical mean-field theory (DMFT), a lattice problem is mapped to an interacting impurity coupled to a self-consistent bath. The exchanged objects are frequency- or time-dependent Green’s functions, hybridization functions, and self-energies. For the standard Hubbard setting, single-site DMFT becomes exact in the infinite-coordination limit with local interactions and hopping scaled to keep the kinetic energy finite (for example, as coordination$^{-1/2}$ for the usual nearest-neighbor constructions), but practical calculations still have impurity-solver, bath-discretization, statistical, and continuation errors.

Density matrix embedding theory (DMET) instead uses a frequency-independent one-particle density matrix and a finite bath derived from an auxiliary state. In standard ground-state DMET, Schmidt decomposition of an auxiliary Slater determinant supplies at most one entangled bath orbital per fragment orbital; this finite construction is not a general exact replacement for a dynamical bath. The closure and observables differ from DMFT. Calling both methods “an impurity loop” is not enough to make their guarantees interchangeable.

A proposed quantum-accelerated embedding workflow can use a quantum processor as the impurity solver while retaining classical density-functional, mean-field, DMFT, or DMET machinery around it. The quantum processor then addresses a selected correlated subproblem. It does not directly simulate all electrons in the bulk material, and the embedding approximation remains part of the scientific error budget.

Fixed-point stability

Write the ideal deterministic update as

$$x_{k+1}=G(x_k),$$

with fixed point $x_\star=G(x_\star)$. Suppose $G$ is a contraction in a declared norm over the relevant region,

$$\left\| G(x)-G(z) \right\| \le L\lVert x-z\rVert, \qquad 0\le L<1,$$

assume both ideal and perturbed iterates stay in that region, and let each implemented update have bounded perturbation $\lVert\xi_k\rVert\le\eta$:

$$\widehat x_{k+1} = G(\widehat x_k)+\xi_k.$$

Then

$$\left\| \widehat x_k-x_\star \right\| \le L^k \left\| \widehat x_0-x_\star \right\| + \eta \frac{1-L^k}{1-L}.$$

The asymptotic error floor is at most $\eta/(1-L)$. Near a weakly stable fixed point with $L\approx1$, small quantum-estimator or fitting errors can be strongly amplified. This bound is illustrative rather than universal: many embedding maps are not globally contractive, may have multiple solutions, and can be non-normal in their linearized dynamics.

Mixing changes stability, not the fixed point

A common damped update is

$$x_{k+1} = (1-\lambda)x_k +\lambda G(x_k), \qquad 0<\lambda\le1.$$

For a scalar map, local stability at $x_\star$ requires

$$\left| 1-\lambda +\lambda G'(x_\star) \right|<1.$$

Smaller $\lambda$ can stabilize an oscillatory update, but it also makes the raw iterate difference

$$\left|x_{k+1}-x_k\right| = \lambda \left|G(x_k)-x_k\right|$$

artificially small. A stopping test should inspect the unmixed physical residual $\lVert G(x_k)-x_k\rVert$, not only successive iterate changes.

In several dimensions, the Jacobian spectrum controls local asymptotic stability, while non-normality can produce large transient amplification even when every eigenvalue lies inside the unit disk. Monitoring only a final fixed-point residual can miss that sensitivity.

Worked example: bias amplification in a scalar closure

Consider

$$G(x)=ax+c, \qquad |a|<1.$$

The ideal fixed point is

$$x_\star = \frac{c}{1-a}.$$

If the quantum solver contributes a persistent additive bias $b$, the implemented map is $\widetilde G(x)=ax+c+b$, and its apparently well-converged fixed point is

$$\widetilde x_\star = x_\star +\frac{b}{1-a}.$$

For $a=0.9$ and $b=10^{-3}$, the final bias is $10^{-2}$. The iteration can converge smoothly to many decimal places while being ten times more biased than the quantum subroutine’s local offset. Convergence certifies consistency with the implemented map, not correctness of that map.

Embedding-specific approximation layers

An embedding result may depend on all of the following:

• choice and size of the correlated fragment or active space;
• bath construction and discretization;
• environment approximation and exchange-correlation functional;
• interaction parameters and double-counting correction;
• impurity state preparation and solver error;
• real- or imaginary-frequency grid and Fourier transforms;
• analytic continuation when real-frequency spectra are inferred;
• statistical error, covariance, and nonlinear fitting;
• branch selection, initialization, mixing, and stopping rule.

Increasing the accuracy of the quantum impurity solver cannot remove errors from the surrounding model. Conversely, a useful embedding approximation may not require a uniformly precise impurity state if the requested fragment observable is robust. Accuracy should be allocated to the final scientific quantity.

Unified Feedback Error Propagation

Variational and embedding workflows can both be written locally as

$$\begin{aligned} \widehat y_k &= Q(x_k)+b_k+\eta_k, \\ x_{k+1} &= F(x_k,\widehat y_k), \end{aligned}$$

where $b_k$ is systematic bias and $\eta_k$ is a random estimator error. Near a reference trajectory, first-order perturbation gives

$$\delta x_{k+1} = J_k\delta x_k +B_k b_k +B_k\eta_k,$$

with

$$J_k = \frac{\partial F}{\partial x} +\frac{\partial F}{\partial y} \frac{\partial Q}{\partial x}, \qquad B_k = \frac{\partial F}{\partial y}.$$

This equation exposes three separate questions:

1. Stability: do perturbations in the current classical state decay or grow under $J_k$?
2. Interface sensitivity: how strongly does the update react to measured quantities through $B_k$?
3. Estimator quality: what are the bias and covariance of the actual quantum output?

For a stationary linearized map $J$ with spectral radius below one, a persistent bias produces the steady displacement

$$\delta x_{\mathrm{bias}} = (I-J)^{-1}B b.$$

For additive zero-mean innovations independent of the current state and of earlier innovations, with finite covariance $\Sigma_\eta$, the stationary state covariance, when it exists, satisfies the discrete Lyapunov equation

$$\Sigma_x = J\Sigma_xJ^{\mathsf T} +B\Sigma_\eta B^{\mathsf T}.$$

The mean and covariance answer different questions. More shots can reduce the random term but not a fixed readout, model, or mitigation bias. Correlated drift violates the independent-noise assumption and can masquerade as physical movement of the classical state.

Digital–analog open-loop sequences do not have the same classical feedback, but coherent block errors propagate through conjugation by subsequent unitaries. In every family, the order and sensitivity of composition matter; subroutine error bars cannot simply be pasted onto the final answer.

Resource Accounting

Hybrid methods often shorten one coherent quantum execution by repeating it many times. A transparent resource statement reports both axes.

For adaptive iteration $k$, let $m_k$ be the number of measurement settings and $n_{ka}$ the shots assigned to setting $a$. The total shots are

$$N_{\mathrm{shots}} = \sum_{k=1}^{K} \sum_{a=1}^{m_k} n_{ka}.$$

This number should be accompanied by:

• quantum systems or logical qubits and any ancillas;
• native analog evolution time and control bandwidth;
• digital gate counts, depth, and two-qubit or multiqubit blocks;
• reset, preparation, readout, and calibration executions;
• number of adaptive round trips and optimizer or embedding restarts;
• classical memory, linear algebra, tensor contractions, fitting, and preprocessing;
• compilation, queue, communication, and data-analysis latency;
• failed, discarded, and postselected executions.

A schematic wall-clock model is

$$T_{\mathrm{wall}} \approx \sum_{k=1}^{K} \Bigg[ T_{\mathrm{compile}}^{(k)} +T_{\mathrm{queue}}^{(k)} +T_{\mathrm{classical}}^{(k)} + \sum_{a=1}^{m_k} n_{ka} T_{\mathrm{cycle}}^{(k,a)} \Bigg].$$

Measurement settings within one iteration may be parallelized, but adaptive iterations generally lie on a critical path. A method with shallow circuits can be slow when it requires thousands of sequential quantum–classical round trips. Conversely, a native analog block may reduce gate count while requiring expensive Hamiltonian characterization over a large control region.

Include the surrounding classical calculation

A claim of quantum benefit should include all classical work that is necessary to construct and interpret the quantum subproblem. For embedding, that may include electronic-structure preprocessing and self-consistency. For variational dynamics, it includes metric estimation, regularization, integration, and repeated compilation. For digital–analog simulation, it includes pulse optimization and generator identification.

The comparison should use the best applicable classical method for the same target observable, accuracy, and input access. A quantum fragment solver must be compared with modern classical impurity solvers, not with brute-force simulation of the entire material. A short trajectory must be compared with tensor-network, Monte Carlo, semiclassical, or problem-specific methods in the regime where they apply.

Validation and Stopping

A hybrid workflow needs tests at three scales:

1. Component tests validate native blocks, circuits, estimators, classical solvers, and update rules independently.
2. Interface tests check basis conventions, units, covariance propagation, serialization, and sensitivity to exchanged quantities.
3. End-to-end tests compare final observables with exact limits, independent methods, held-out data, or experimentally known constraints.

Validation interventions

Useful interventions include:

• solve small instances exactly and compare the entire trajectory, not only its endpoint;
• halve a time step, tighten a regularizer, enlarge an ansatz, or increase a fragment and look for controlled convergence;
• replace measured quantum data with exact simulated interface data to isolate the classical loop;
• replay recorded quantum data through alternative classical updates;
• test symmetries, conservation laws, positivity, causality, and sum rules;
• reserve observables that are not used for fitting or stopping;
• interleave reference circuits or calibration blocks to expose drift;
• compare branches from multiple initializations and report the selection rule;
• perturb interface data within its covariance and propagate the full result.

For digital–analog protocols, turning off selected refocusing pulses can test the inferred Hamiltonian terms. For variational dynamics, compare the measured projection residual with held-out observables and step-size convergence. For embedding, vary fragment size, bath representation, mixing, and initial branch while monitoring physical residuals.

A stopping rule is not a certificate

A robust stopping rule combines several conditions, for example:

$$\begin{aligned} \lVert r_k\rVert &\le \tau_r, \\ \left| \widehat O_k-\widehat O_{k-1} \right| &\le \tau_O, \\ \operatorname{SE}(\widehat O_k) &\le \tau_{\mathrm{stat}}, \end{aligned}$$

for a required number of consecutive iterations, together with a maximum budget. The residual $r_k$ must be defined in physical, unmixed variables. Repeated estimates in an adaptive run can be correlated through their settings or shared calibration. The standard error of a difference must include any such covariance. Repeatedly checking ordinary fixed-sample confidence intervals until a favorable stopping event can also invalidate their nominal coverage; a sequential guarantee or an independent final assessment is needed for that claim.

Meeting these conditions establishes that the implemented loop has stopped according to its declared rule. It does not establish low model error, uniqueness of the solution, or agreement with the target system. Those require validation evidence beyond the stopping data.

Choosing a Hybrid Architecture

Scientific or engineering constraint	Architecture that may help	Main price paid
a calibrated many-body interaction is much cheaper than decomposing it	digital–analog blocks	hardware-specific synthesis and generator uncertainty
target dynamics remain near a compact parameterized manifold	variational projection	repeated metric and force estimation, projection bias, conditioning
strong correlations are localized to a fragment or active space	quantum embedding	environment approximation, boundary matching, self-consistency
coherence time is scarce but repeated preparation is available	variational or embedding loop	shot count, latency, and feedback bias
model parameters must change widely	digital gates or variational controls	compilation or retraining cost
direct access to large analog systems is available	analog or digital–analog protocol	restricted model family and validation burden

The architecture should follow the target bottleneck. Hybridization for its own sake adds interfaces, and every interface is another place where conventions, uncertainty, or causality can be lost.

A useful architecture addresses a particular bottleneck; hybridization itself supplies no accuracy or complexity theorem. Digital Quantum Simulation research directions discusses open problems in long-time variational accuracy and experimentally reliable simulation methods.

Exercises

1. Classify interfaces, not labels

Classify each workflow and identify the quantity crossing its principal interface:

1. fixed single-qubit rotations alternate with tunable all-to-all Ising evolution;
2. a circuit estimates a tangent-space metric and a classical integrator updates circuit angles;
3. a quantum processor returns an impurity Green’s function used to update a classical bath;
4. a classically optimized pulse implements a fixed gate and is then used without feedback in an otherwise digital circuit.

Solution

1. This is digital–analog simulation. The interface inside the quantum control sequence exchanges frame choices, pulse labels, and native-block durations.
2. This is variational quantum–classical simulation. Estimated entries of $A$ and $\boldsymbol C$, with covariance, cross from quantum execution to the classical integrator; updated angles return.
3. This is quantum embedding. A Green’s function crosses to the classical bath update, and new bath or hybridization parameters return.
4. The offline use of classical pulse optimization does not by itself make the simulation hybrid in the architectural sense used here. Every implemented gate has classical control. The executed simulation is digital unless the optimized pulse is exposed as a continuously parameterized analog block or enters an adaptive quantum–classical loop.

2. Verify the two-qubit echo

For

$$H=JZ_1Z_2+h_1Z_1+h_2Z_2,$$

show explicitly that the echo sequence in the text cancels both local fields. Would $P=X_1$ cancel the same terms?

Solution

Conjugation by $P=X_1X_2$ gives

$$PZ_1P=-Z_1, \qquad PZ_2P=-Z_2, \qquad PZ_1Z_2P=Z_1Z_2.$$

Because all three Pauli strings commute,

$$U_{\mathrm{echo}} = \exp\!\left[ -\frac{i\tau}{2\hbar}(PHP+H) \right] = e^{-iJZ_1Z_2\tau/\hbar}.$$

For $P=X_1$, both $Z_1$ and $Z_1Z_2$ change sign while $Z_2$ does not. That echo would cancel the interaction and the $h_1Z_1$ term but retain $h_2Z_2$.

3. Derive the leading composition term

Use the Baker–Campbell–Hausdorff formula to find the leading correction to the effective Hamiltonian of $e^{-iA\Delta t/\hbar}e^{-iB\Delta t/\hbar}$. Explain why the correction is Hermitian when $A$ and $B$ are Hermitian.

Solution

Set $X=-iA\Delta t/\hbar$ and $Y=-iB\Delta t/\hbar$. Then

$$\log(e^Xe^Y) = X+Y+\frac12[X,Y]+O(\Delta t^3).$$

Since

$$[X,Y] = -\frac{\Delta t^2}{\hbar^2}[A,B],$$

comparison with $-iH_{\mathrm{eff}}\Delta t/\hbar$ gives

$$H_{\mathrm{eff}} = A+B -\frac{i\Delta t}{2\hbar}[A,B] +O(\Delta t^2).$$

For Hermitian $A$ and $B$, the commutator is anti-Hermitian, so $-i[A,B]$ is Hermitian.

4. Derive the projected real-time equation

Starting from

$$\lvert R\rangle = \sum_j\dot\theta_j\lvert D_j\psi\rangle +\frac{i}{\hbar}(H-E)\lvert\psi\rangle,$$

minimize $\langle R\vert R\rangle$ with respect to real $\dot\theta_i$. Recover $A\dot{\boldsymbol\theta}=\boldsymbol C$.

Solution

For real velocities,

$$\frac12 \frac{\partial}{\partial\dot\theta_i} \langle R\vert R\rangle = \operatorname{Re} \langle D_i\psi\vert R\rangle.$$

Stationarity gives

$$\sum_j \operatorname{Re} \langle D_i\psi\vert D_j\psi\rangle \dot\theta_j +\operatorname{Re} \left[ \frac{i}{\hbar} \langle D_i\psi\rvert H\lvert\psi\rangle \right] =0.$$

The energy term vanishes because $\langle D_i\psi\vert\psi\rangle=0$. Using $\operatorname{Re}(iz)=-\operatorname{Im}(z)$ gives

$$\sum_j A_{ij}\dot\theta_j = \frac{1}{\hbar} \operatorname{Im} \langle D_i\psi\rvert H\lvert\psi\rangle =C_i.$$

5. Diagnose an ill-conditioned variational metric

Let

$$A= \begin{pmatrix} 1&0\\ 0&10^{-6} \end{pmatrix}, \qquad \boldsymbol C= \begin{pmatrix} 1\\ 10^{-3} \end{pmatrix}.$$

Find the unregularized velocity. If the uncertainty in $C_2$ is $2\times10^{-4}$, find its contribution to the uncertainty in $\dot\theta_2$. What changes if eigenvalues below $10^{-5}$ are truncated?

Solution

Direct inversion gives

$$\dot{\boldsymbol\theta} = \begin{pmatrix} 1\\ 10^3 \end{pmatrix}.$$

The standard uncertainty propagated from $C_2$ alone is

$$\sigma_{\dot\theta_2} = \frac{2\times10^{-4}}{10^{-6}} =200.$$

The small metric eigenvalue amplifies both the force and its noise. Truncating eigenvalues below $10^{-5}$ sets the second pseudoinverse component to zero, so the retained solution is $(1,0)^{\mathsf T}$. Its variance is lower, but it has deliberately removed motion along the second direction. Whether that bias is acceptable must be tested against trajectory observables.

6. Prove the contraction error floor

Assume $G$ has contraction factor $L<1$ and

$$\widehat x_{k+1}=G(\widehat x_k)+\xi_k, \qquad \lVert\xi_k\rVert\le\eta.$$

Derive the finite-$k$ error bound relative to the ideal fixed point.

Solution

Let $e_k=\lVert\widehat x_k-x_\star\rVert$. Since $x_\star=G(x_\star)$,

$$e_{k+1} \le L e_k+\eta.$$

Iterating the recursion gives

$$e_k \le L^k e_0 +\eta\sum_{j=0}^{k-1}L^j = L^k e_0 +\eta\frac{1-L^k}{1-L}.$$

Thus $\limsup_{k\to\infty}e_k\le\eta/(1-L)$. The result requires a uniform deterministic perturbation bound, contraction in the chosen norm, and iterates remaining in its domain; it should not be applied unchanged to unstable or multibranch maps.

7. Analyze damping and a deceptive residual

For $G(x)=-2x+c$, determine the range of $\lambda$ for which

$$x_{k+1}=(1-\lambda)x_k+\lambda G(x_k)$$

is locally stable. Explain why $|x_{k+1}-x_k|$ alone becomes a poor stopping diagnostic as $\lambda\to0$.

Solution

The derivative of the mixed map is

$$1-\lambda+\lambda(-2) = 1-3\lambda.$$

Stability requires

$$|1-3\lambda|<1,$$

so $0<\lambda<2/3$. Meanwhile,

$$x_{k+1}-x_k = \lambda\bigl(G(x_k)-x_k\bigr).$$

For very small $\lambda$, successive iterates can be nearly equal even when the unmixed fixed-point residual $G(x_k)-x_k$ is large. The stopping rule should examine the latter or divide out the known mixing factor.

8. Correlated measurement drift in a stable loop

Consider $x_{k+1}=a x_k+\eta_k$, with $|a|<1$. The stationary measurement error has variance $\sigma^2$ and correlations $\mathbb E[\eta_k\eta_{k-j}]=\sigma^2\rho^{|j|}$, $|\rho|<1$. Derive the stationary variance of $x$. Compare $a=0.9,\rho=0.8$ with independent noise of the same variance.

Solution

The stationary solution is $x_k=\sum_{j\geq0}a^j\eta_{k-1-j}$, so

$$\mathbb E[x_k\eta_k] =\frac{\rho\sigma^2}{1-a\rho}.$$

Taking the variance of the recursion and retaining this cross term gives

$$\operatorname{Var}(x) =\frac{\sigma^2}{1-a^2} \frac{1+a\rho}{1-a\rho}.$$

Independent noise gives $\sigma^2/(1-a^2)$. At the stated values, persistent drift increases this by $1.72/0.28\simeq6.14$. Both processes have the same single-step variance; their different time correlations produce different final uncertainties. Extra repeated shots within a shared drift interval need not remove that drift.

References

• I. Arrazola, J. S. Pedernales, L. Lamata, and E. Solano, “Digital-analog quantum simulation of spin models in trapped ions,” Scientific Reports 6, 30534 (2016), doi:10.1038/srep30534.
• B. Bauer, D. Wecker, A. J. Millis, M. B. Hastings, and M. Troyer, “Hybrid quantum-classical approach to correlated materials,” Physical Review X 6, 031045 (2016), doi:10.1103/PhysRevX.6.031045.
• S. Endo, J. Sun, Y. Li, S. C. Benjamin, and X. Yuan, “Variational quantum simulation of general processes,” Physical Review Letters 125, 010501 (2020), doi:10.1103/PhysRevLett.125.010501.
• A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, “Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions,” Reviews of Modern Physics 68, 13–125 (1996), doi:10.1103/RevModPhys.68.13.
• G. Knizia and G. K.-L. Chan, “Density matrix embedding: a simple alternative to dynamical mean-field theory,” Physical Review Letters 109, 186404 (2012), doi:10.1103/PhysRevLett.109.186404.
• Y. Li and S. C. Benjamin, “Efficient variational quantum simulator incorporating active error minimization,” Physical Review X 7, 021050 (2017), doi:10.1103/PhysRevX.7.021050.
• S. McArdle, T. Jones, S. Endo, Y. Li, S. C. Benjamin, and X. Yuan, “Variational ansatz-based quantum simulation of imaginary time evolution,” npj Quantum Information 5, 75 (2019), doi:10.1038/s41534-019-0187-2.
• A. D. McLachlan, “A variational solution of the time-dependent Schrödinger equation,” Molecular Physics 8, 39–44 (1964), doi:10.1080/00268976400100041.
• A. Parra-Rodriguez, P. Lougovski, L. Lamata, E. Solano, and M. Sanz, “Digital-analog quantum computation,” Physical Review A 101, 022305 (2020), doi:10.1103/PhysRevA.101.022305.
• X. Yuan, S. Endo, Q. Zhao, Y. Li, and S. C. Benjamin, “Theory of variational quantum simulation,” Quantum 3, 191 (2019), doi:10.22331/q-2019-10-07-191.