Digital quantum simulation of open quantum systems using quantum imaginary time evolution

I. INTRODUCTION

• Coupling to the environment leads to non-unitary evolution, which is not naturally simulable on quantum hardware; QITE, introduced for ground-state preparation and thermal averages on near-term devices, is adapted for open-system evolution.Source: §I · pp. 2–3
• Existing approaches can require known Kraus operators, system-environment evolution, or operator-sum representations, while variational methods require an ansatz and potentially high-dimensional classical optimization.Exactly determining Kraus operators from a Lindblad equation is classically hard and can be applied only to small systems; explicit Trotterization has ancilla overhead, while variational optimization is NP-hard.Source: §I · pp. 2–3
• The first algorithm vectorizes the density operator, Trotterizes its Hermitian evolution, and uses QITE for the anti-Hermitian component.The Lindblad equation is cast as a Schrödinger-type equation with a non-Hermitian Hamiltonian.Source: §I · p. 3
• The second algorithm uses an ansatz for the density operator that is preserved during real- and imaginary-time evolution.Source: §I · p. 3
• The hardware results show good agreement with exact results for spontaneous emission and the dissipative transverse field Ising model.The examples are a two-level system in a zero-temperature heat bath and a two-site model.Source: §I · p. 3

II. THEORY

• The algorithms target Lindblad-form dynamics, whose derivation is often based on weak system-environment coupling and absence of memory effects.These are the Born-Markov approximation assumptions.Source: §II · pp. 3–4
• Figure 1 depicts Algorithm I as vectorized-density evolution and Algorithm II as purification-based evolution, with non-unitary terms applied to the system for a time step.The figure text also identifies a generalized 2n-qubit Bell-state preparation in Algorithm I and a qubit bundle of size n.Source: Figure 1 · p. 4, §II · p. 4
• Algorithm I vectorizes the density operator and Trotterizes its evolution, while Algorithm II combines a QITE adaptation with an ansatz for the time-dependent density operator.Source: §II · p. 4

A. Algorithm I

• Algorithm I maps the density operator to a vector in an extended Hilbert space and expresses Lindblad evolution as Schrödinger-type evolution with a non-Hermitian Hamiltonian.The density operator is column-stacked before the transformation.Source: A. Algorithm I · pp. 4–5
• The non-unitary propagator exp(−H2τ) is implemented by analogy to QITE, while the Hermitian component is applied as unitary evolution.The evolution is Trotterized with time step τ = t/N.Source: A. Algorithm I · p. 5
• The algorithm preserves Tr(ρ²), not Tr(ρ), so the matricized output is not strictly a density matrix and observables require trace normalization.The final expectation value is normalized as ⟨O⟩/Tr(ρ), requiring both quantities at each time step.Source: A. Algorithm I · p. 6
• Observable measurements use a Hadamard test to obtain unnormalized expectation values, with trace evaluation also required for normalized results.For k-local observables, the measurement construction uses a generalized Bell state on n−k qubits.Source: A. Algorithm I · pp. 6–7

B. Algorithm II

• The density-operator ansatz sums weighted computational-basis states over an index set I, with a strict subset giving a sparse approximation.The representation is exact only when I is large enough; polynomial resources for storing weights imply sparsity.Source: B. Algorithm II · pp. 7–8
• The Lindblad propagator is Trotterized term by term, and its non-unitary terms are treated using QITE-style approximations involving weights qx and a Hermitian operator A.The text describes separate forms for the two non-unitary terms.Source: B. Algorithm II · pp. 8–9
• Algorithm II requires no ancilla qubits or Hadamard test for observable measurements, and can halve the required qubit count while trading for more classical and quantum computation time.The text identifies this resource trade-off as useful for near-term hardware.Source: B. Algorithm II · p. 9
• The authors identify the sparse representation of the density matrix and the number of measurements required to evolve the system as drawbacks of the approach.Source: B. Algorithm II · p. 9

III. RUN TIME BOUNDS, COMPUTATIONAL OVERHEADS, AND ERRORS

• Hardware-noise errors are not addressed because the analysis treats them as non-algorithmic.The section frames its analysis around runtimes, computational overheads, and algorithm-associated errors.Source: §III · pp. 9–10

A. Run time bounds

• The runtime analysis assumes a local Lindblad equation and treats the imaginary-time propagator as an additional Trotterized evolution.The propagator exp(−H2τ) is divided into p2 steps for a desired error ϵ2.Source: A. Run time bounds · p. 10
• For the imaginary-time approximation, p2 > 1/ϵ2 bounds the error by ϵ2 when 3m2tv2/N < 1.Here v2 = max_l{∥hl∥}.Source: A. Run time bounds · p. 10
• Algorithm II has the same runtime as Algorithm I up to constant factors, but its bound excludes errors from approximating the density operator with a strict subset of bit strings.Excluding bit strings creates representation error that would need assessment by stochastic sampling.Source: A. Run time bounds · pp. 10–11

B. Measurement and classical computational overheads

• Both algorithms incur measurement overhead in addition to running time, and their asymptotic per-step circuit requirements are summarized for open systems on n sites.The table gives counts per time step per Lindblad operator.Source: B. Measurement and classical… · p. 11, Table I · p. 12
• Algorithm I retains the original QITE measurement overhead: state tomography over each D-qubit domain constructs the unitaries, with O(n/D) domains in a one-dimensional lattice.The supplied passage also gives a measurement count fragment, but it does not support a reliable complete expression.Source: B. Measurement and classical… · p. 11
• For efficiency, the subset I of bit strings whose matrix elements are measured must scale polynomially or slower; measuring all matrix elements requires runs of circuits.The asymptotic circuit count also depends on the number of Lindblad operators and the size of I.Source: B. Measurement and classical… · p. 11
• Finite-domain approximation is accurate in many cases, and dissipation can permit a smaller domain D than in the same system without dissipation because it generally reduces correlation length.The cost-reducing correlation length may also permit an efficient classical description of the evolution.Source: B. Measurement and classical… · p. 11
• Whether reduced correlation length also enables efficient classical description marks an imprecise boundary for quantum simulation algorithms and remains under active investigation.The passages identify this as a general consideration, not a settled threshold.Source: B. Measurement and classical… · p. 11

IV. RESULTS

• For TLS observables, hardware results from both algorithms qualitatively agree with theory, with deviations largely attributable to gate errors.The experiments used IBM hardware; the population curves approach a non-equilibrium steady state for γt ≳ 5.Source: §IV · pp. 13–14, Figure 2 · p. 13
• The TLS exhibits initial exponential excited-state decay, damped Rabi oscillations, and a non-equilibrium steady state where total emission balances absorption.The steady-state populations reflect bath-induced spontaneous emission balanced against absorption and stimulated emission from Hamiltonian driving.Source: §IV · p. 14
• At the driven TLS steady state, off-diagonal coherences remain nonzero, while normalized purity varies over time despite preservation of the simulator’s inner product.The persistent coherence indicates the hardware captures quantum coherence in the canonical basis.Source: §IV · p. 14
• For the two-site dissipative TFIM, both algorithms qualitatively match exact magnetization dynamics; Algorithm I agrees at Trotter step γt/N ∼ 0.5, and Algorithm II tracks all simulated times.Magnetization oscillations arise from the relatively large transverse field.Source: Figure 3 · p. 15, §IV · pp. 15–16
• Algorithm I accurately describes the dissipative dynamics using 16 Pauli strings; Fig. S2 shows no significant increase in accuracy when using more than 16 strings.The required number of Pauli strings is problem-dependent.Source: §IV · p. 16
• Each TFIM time step requires 836 circuits for Algorithm II to measure two-qubit Pauli-string matrix elements, versus 16 circuits for Algorithm I’s 16 operators.Algorithm II requires measuring all two-qubit Pauli-string matrix elements at each time step.Source: §IV · p. 16

• At larger dissipation rates γ ∼ J, h, both algorithms accurately capture magnetization dynamics in separate simulations (Fig. S3).Source: §IV · p. 16
• The simulations omit Algorithm II’s error from using a subset of bit-strings, so its actual error depends on how accurately that subset represents the density matrix for the problem.Source: §IV · p. 16
• Stochastic sampling of bit-strings may be viable for larger systems.Source: §IV · p. 16

V. SUMMARY

• The paper introduces digital quantum algorithms based on QITE for time evolution of open quantum systems governed by a Lindblad equation.Source: §V · pp. 16–17
• Algorithm I applies QITE to non-unitary evolution after vectorizing the density operator.Source: §V · p. 16
• Algorithm II adapts QITE to preserve a purification-based ansatz during computation.Source: §V · p. 16
• IBM Quantum processor simulations showed good qualitative agreement with exact results in all cases and reduced the resources required for Lindblad-system simulations.The demonstrations covered spontaneous emission in a two-level system and the dissipative transverse-field Ising model.Source: §V · pp. 16–17

Supplemental Materials: Digital quantum simulation of open

• The supplemental-material title identifies quantum imaginary time evolution as the approach for digital quantum simulation of quantum systems.Source: Supplemental Materials: Digital… · p. 23
• The listed authors are Hirsh Kamakari, Shi-Ning Sun, Mario Motta, and Austin J. Minnich.Source: Supplemental Materials: Digital… · p. 23
• The affiliations listed are Caltech and IBM Quantum, and the document is dated November 18, 2021.Source: Supplemental Materials: Digital… · p. 23

A. Implementing Vk ⊗Vk via a QITE adaptation

• For V_k ⊗ V_k, Algorithm II finds real coefficients q_x and a Hermitian operator A by solving the linear system S a = b.A is decomposed into Pauli strings acting on at most D qubits, with coefficients a_j = O(τ); the optimal coefficients minimize the stated objective.Source: A. Implementing Vk ⊗Vk via a QITE… · pp. 24–25

B. Implementing Wk via a QITE adaptation

• For W_k = exp(τL_k ⊗ L_k), the QITE equations are analogous to those for V_k ⊗ V_k, with the same S matrix but a modified b vector.The vectors |v_j⟩ are identical in both cases, while |Φ⟩ differs.Source: B. Implementing Wk via a QITE adaptation · p. 26

C. Conservation of Probability

• Algorithm II preserves the trace when the q_x contributions from both V_k and W_k sum to zero at each time step.Source: C. Conservation of Probability · p. 26

D. Measuring Observables

• The evolved weights p_x(t) and Hermitian operator A(t) yield observable expectations, but storing every p_x becomes infeasible beyond a certain system size.The text proposes stochastic sampling; locality conditions suggest it can converge faster than uniform random sampling.Source: D. Measuring Observables · pp. 26–27

E. Measuring Matrix Elements

• Matrix elements needed for q_x and a_j can be measured by decomposing operators into Pauli strings and preparing states based on the bitstrings x and y.For a k-qubit local observable, only the first k qubits need a superposition-state preparation, which can be efficient when k is small; locality-guided sampling favors the first k + m qubits, and dissipation may reduce long-range correlations and improve convergence.Source: E. Measuring Matrix Elements · pp. 27–28

II. HARDWARE ERRORS IN THE TWO LEVEL SYSTEM SIMULATION

• Algorithm II's hardware results deviated less from noiseless emulation than Algorithm I's, attributed to Algorithm I's greater circuit depth and qubit count.Algorithm II's two-level-system simulation required only a single qubit and single-qubit rotations, avoiding infidelity from two-qubit gates.Source: Figure S1 · p. 29, §II · p. 29

III. NUMBER OF PAULI STRINGS IN THE TRANSVERSE FIELD ISING

• Including more than 16 Pauli strings produced only a minor effect on the transverse-field Ising model observables.The study compared simulations with increasing Pauli-string counts against the exact result at dissipation rate γ = 0.1.Source: §III · p. 30, Figure S2 · p. 30

IV. EFFECT OF DISSIPATION RATE ON ALGORITHM PRECISION

• Larger dissipation rates correspond to larger Lindblad-operator norms, which increase both Trotter and QITE errors.Source: §IV · p. 31
• Algorithm II qualitatively outperformed Algorithm I at every simulated dissipation rate in the 2-site TFIM, using all bit-strings versus 16 Pauli strings.Source: §IV · p. 31, Figure S3 · p. 31
• The simulations did not account for bit-string-selection error, which would add error in larger systems where all bit-strings cannot be included.Source: §IV · p. 31