Adaptive variational simulation for open quantum systems

Emerging quantum hardware provides new possibilities for quantum simulation. While much of the research has focused on simulating closed quantum systems, the real-world quantum systems are mostly open. Therefore, it is essential to develop quantum algorithms that can effectively simulate open quantum systems. Here we present an adaptive variational quantum algorithm for simulating open quantum system dynamics described by the Lindblad equation. The algorithm is designed to build resource-efficient ansatze through the dynamical addition of operators by maintaining the simulation accuracy. We validate the effectiveness of our algorithm on both noiseless simulators and IBM quantum processors and observe good quantitative and qualitative agreement with the exact solution. We also investigate the scaling of the required resources with system size and accuracy and find polynomial behavior. Our results demonstrate that near-future quantum processors are capable of simulating open quantum systems.


Introduction
The theory of open quantum system investigates the behavior of small quantum systems in contact with a large environment [1][2][3].Recently there has been a surge of interest in developing efficient simulation algorithms for open quantum systems [4][5][6][7][8].On one hand, the emergence of quantum technologies for building artifi-Huo Chen: huochen@lbl.govNiladri Gomes: niladri@lbl.gov,H.C. and N.G.contributed equally.
Siyuan Niu: siyuanniu@lbl.govWibe Albert de Jong: WAdeJong@lbl.govcial quantum systems for computing [9,10], sensing [11,12], light harvesting [13,14], and other applications promises an abundant scientific and societal benefits.These artificial quantum systems must be treated and designed as open because perfect isolation of a quantum system is extremely challenging, if not impossible.On the other hand, the theory of open quantum systems plays a crucial role in the fundamental science because it describes many non-equilibrium processes, including quarkonium suppression in heavy-ion collisions [15], relaxation dynamics of many-body quantum system [4], charge and energy transfer dynamics in molecular systems [16], and others.All of these applications will benefit from improved simulation algorithms for open quantum systems.Similar to the case of closed system, the complexity of classical algorithms for simulating open system scales exponentially with the system size.As the quantum computers hold the promise of being able to efficiently simulate quantum systems, it is worth exploring the new opportunities offered by the rapid advancement of cloudaccessible quantum computers [17][18][19][20][21][22][23][24][25][26] for opensystem simulations.Current research mostly is focused on solving the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) equation (or simply Lindblad equation) due to its broad applications in various subfields of quantum systems engineering [1,3,[27][28][29].Based on the underlying techniques, existing quantum algorithms for simulating the Lindblad dynamics fall into five categories: unitary dilation [6,18,30,31], variational simulation [5], quantum imaginary-time evolution (QITE) [32], Monte Carlo [33] and qubitbath engineering [19,34,35].Among them the variational simulation and QITE require less circuit depth and are considered more NISQ (noisy intermediate-scale quantum) [36] friendly.
Variational algorithms have been shown to be capable of solving real and imaginary time evolutions [37][38][39], excited states [40], thermal states [41], non-Hermitian quantum mechanics for nonequilibrium systems [42], open quantum systems and general first-order differential equations [5].In this paper, we present a compact approach for solving the Lindblad equation using a timedependent adaptive variational method.Our strategy to generate NISQ-friendly ansätze is built upon adaptive variational quantum dynamics simulations (AVQDS) [43] and its imaginarytime counterpart (AVQITE) [39].Distinct from previous works, our focus is on the dynamics of open quantum systems, rather than the closed system dynamics or ground state solutions typically explored.The main idea is to reformulate the Lindblad equation as a Schrödinger equation with an effective non-Hermitian Hamiltonian, and utilizing a quantum computer to simulate the evolution of the normalized state vector while employing a classical computer to record its norm.The presence of both unitary and dissipative elements in the equation complicates the adaptive process, as it lacks a pre-defined lower bound, rendering the establishment of a practical threshold impractical.To address this challenge, we have developed an unrestricted adaptive protocol, specifically tailored to circumvent the difficulties in setting a feasible threshold for the McLachlan distance.
To showcase the efficacy of our algorithm, we simulate the open quantum system dynamics of a quantum annealing (QA) process, focusing on the alternating-sector-chain (ASC) problem as described in [44,45].QA is a quantum meta-heuristic algorithm that is commonly implemented on analog quantum hardware.Although our algorithm is specifically designed for gatebased quantum computers, simulating QA hardware highlights its applicability to a wider range of use cases, such as analyzing artificial quantum systems.Furthermore, the ASC problem stands out as an effective toy problem to benchmark open-system algorithms.On one hand, the system is integrable in the closed-system case because its Hamiltonian can be transformed into a free-fermion model by Jordan-Wigner transformation.On the other hand, the introduction of local dephasing and amplitude damping noise perturb the free-fermion model and make it a non-trivial task to simulate the quantum dynam-ics.Given that dephasing and amplitude damping represent two of the most prevalent quantum noise models, holding significant relevance in realworld applications [30,32,33], the open-system ASC problem serves as a ideal test platform for open-system algorithms.
We demonstrate our algorithm on both the noiseless simulator and IBM quantum hardware.Our results show good quantitative and qualitative agreement with the exact solution.Furthermore, we perform numerical analysis of the resource requirement to simulate small systems, and find polynomial scaling with respect to either the system size or the desired accuracy.These results provide compelling evidence that open quantum system simulation is within the capability of near-term quantum devices.
In summary, we propose a new adaptive variational quantum algorithm tailored for simulating the Lindblad master equation, where prior adaptive approaches encounter the challenge of establishing a feasible threshold for the McLachlan distance due to the coexistence of unitary and dissipative components.We validate the effectiveness of our algorithm on both simulators and IBM quantum hardware.Our results show that open quantum system simulation is within the capability of near-term quantum devices.The structure of this paper is as follows.In Sec. 2, we introduce the Lindblad equation and describe our quantum algorithm, unrestricted adaptive variational quantum dynamics (UAVQDS), for solving it.In Sec.3.1 we report the performance and resource requirement of our algorithm based on the results from the noiseless simulator.In Sec.3.2, we report the performance of our algorithm on IBM quantum computers.We conclude in Sec. 4, and present additional technical details and derivations in the Appendices.

Method
For ease of reference, a table summarizing frequently used symbols can be found in Appendix A.

The Lindblad master equation
When a quantum system interacts with its environment, we can use the Lindblad master equation (ME) to model its behavior.Unlike closed quantum systems where a single quantum state (pure state) can describe the system's complete information, open quantum systems require a density matrix, which is a statistical ensemble of multiple quantum states (mixed state) to fully describe the system's behavior.Here, we focus on the Lindblad equation in the diagonal form where ρ(t) is the density matrix and H is the Hamiltonian describing the system of interest.
The environment density matrix is traced out while deriving the Lindblad equation and its interaction with the system (in the weak coupling regime) is taken care of by a dissipative term L[ρ(t)] which is given by the form, (2) The operators L k are known as the Lindblad operators, and are taken to be dimensionless.γ k are the dissipation rates, which are non-negative and possess a dimension of inverse time.Usually it is convenient to absorb γ k into the definition of Throughout our discussions, we have adopted the convention of setting ℏ = 1.
Because a quantum computer only evolves a pure state under unitaries, the density matrix equation (Eq.(1)) cannot be directly solved by a gate-based quantum computer.This challenge can be addressed by two different approaches: 1) vectorizing [46] or 2) stochastic unravelling of the Lindblad ME [1,[47][48][49].We briefly review here how these two methods work.
Vectorization -The core idea of vectorization is to convert the density matrix to a vector by stacking the columns of ρ on top of one another: where D × D is the dimension of ρ.We will also use the symbol |ρ⟫ interchangeably with vec(ρ).
The map vec is a linear isometry between the D × D Liouville space (Hilbert space under the Hilbert-Schmidt inner product) of ρ and the D 2 Hilbert space of |ρ⟫, with the preserved norms being the trace norm and L 2 norm where H vec is an effective non-Hermitian Hamiltonian given by (see Appendix B for the derivation) where Lk denotes the conjugate of L k and I is the identity matrix matching the dimension of H.We will refer to this approach as the vectorization method for the remainder of the paper.The vectorization method is a well-established technique for the classical simulation of open quantum systems.
However, its application within quantum algorithms remains largely underexplored.To the best of our knowledge, reference [32] stands as the sole work discussing this methodology in the context of simulating open system dynamics on quantum computers.In this cited work, the authors employ QITE to solve the vectorization Lindblad equation and chose ansätze based on numerical experiments.Contrarily, our study presents a systematic approach to construct resource-efficient ansätze, enabling us to execute larger simulations on IBM quantum computers.
An alternative way to vectorize the ME is to expressing the density matrix in the basis of Pauli matrices and representing the superoperators through the Pauli transfer matrices (PTMs).However, a straightforward application of this technique could significantly increase the circuit complexity, making it less suitable for demonstrations on NISQ hardware.We provide an discussion of this approach in Appendix C.
Stochastic unravelling -The Lindblad master equation can be unravelled using quantum trajectory method [1,[47][48][49].For each trajectory, the evolution is governed by a deterministic evolution and a jump process.The deterministic evolution is described by the Schrödinger equation associated with a non-Hermitian Hamiltonian where ψ(t) is the unnormalized state vector and the effective Hamiltonian is given by For an infinitesimal time step from t to t + dt, there are two possible evolutions for ψ(t) : either the deterministic evolution subject to Eq. (8) (with probability 1 − dp) or a jump occurring (with probability dp).The probability of a jump is given by (10) Furthermore, if a jump happens, the unnormalized state is updated as (12) Let us denote the normalized state vector by |ψ j (t)⟩ = ψj (t) / ψj (t) ψj (t) where the index j represent the jth trajectory.The density matrix solution to Eq. (1) can be constructed by ρ(t) = 1 n n j=1 |ψ j (t)⟩⟨ψ j (t)| for large enough n.We will refer to this approach as the trajectory method henceforth.

Solving the effective Schrödinger equation
The key step in both the aforementioned methods is to solve an effective Schrödinger equation with a non-Hermitian Hamiltonian, which can be accomplished using variational quantum algorithms [50,51].Without loss of generality, we use H eff to represent the effective Hamiltonian in both methods, and ψ(t) , |ψ(t)⟩ to represent the unnormalized and normalized state vector solutions.Furthermore, we split H eff into its Hermitian and anti-Hermitian parts and denotes them by H e and H a respectively Inspired by AVQDS [43], we propose a similar adaptive protocol for the open quantum system simulation.We briefly describe the non-adaptive algorithm here and leave the adaptive protocol to next section.The core idea of an variational algorithm is to encode quantum states with a sequence of parameterized circuits where e −iθµ(t)Aµ is the µth layer of circuit controlled by the real parameter θ µ (t), and |ψ R ⟩ is a fixed reference state.Then the evolution of the state can be mapped to the evolution of the parameters controlling the circuit, i.e., θ µ (t) in Eq. (14).For adaptive protocols, the ansatz operators A µ are adaptively added from an ansatz pool at runtime.Our algorithm is based on McLachlan variational principle [52] where The norm in Eq. ( 15) is known as the McLachlan distance and we denote its square by D. The corresponding evolution equation of the variational parameters is given by where M is a matrix with elements ) and V is a vector with elements where ⟨H eff ⟩ = ⟨ϕ|H eff |ϕ⟩.Detailed derivation of the evolution equation can be found in the Appendix D. Assuming every ansatz operator is a single Pauli string, each term in Eqs.(17) and (18) involving the derivatives can be evaluated using a quantum computer by either direct or indirect measurements [43,53] (See Appendix E for a description of the relevant quantum circuits).
Before we proceed, it is important to note that our algorithm effectively utilizes only the ansatz state (Eq.( 14)) to track the evolution of the normalized state vector.We record the evolution of the state vector norm using classical memory.As detailed in Appendix F, for an infinitesimal time step dt, the norm of ψ(t) shrinks according to where dΓ = 2 ⟨ψ(t)|H a |ψ(t)⟩ dt, which is equivalent to ⟨ϕ|H a |ϕ⟩, can be measured using the same circuit that is used to measure the effective Hamiltonian.In Fig. 1, we provided a single qubit example using Bloch sphere to demonstrate of how our algorithm works.

Restricted and unrestricted adaptive protocol
Unlike existing variational methods [5] that use a fixed set of ansatz operators, our approach uses an adaptive procedure that selects ansatz operators from a predefined pool during the time evolution.The purpose of the adaptive protocol is to define a systematic way of preparing a compact ansatz.This is achieved by minimizing the McLachlan distance [52] during each time step such that there are always enough operators in the ansatz to keep the distance below a threshold.However, a unique challenge for adaptive protocols is the lack of a known attainable lower bound for the McLachlan distance a priori, making it difficult to set a threshold value.This is an issue for both closed and open system adaptive algorithms.For the closed system case, if there are infinite operators in the pool, the lower bound should, in principle, be 0.However, the same is not true for the open system case, which we will explain later.To address this issue, we propose the unrestricted adaptive variational quantum dynamics (UAVQDS) protocol, where instead of setting a fixed threshold for the McLachlan distance, we adopt a greedy approach which selects every operator that lowers the McLachlan distance.
We formulate two versions of the adaptive protocol, the restricted version and unrestricted version.In the restricted version, the McLachlan's distance is set at a fixed threshold.While evolving, at every time step we measure the McLachlan's distance and a number bigger than the threshold triggers the adaptive procedure.We will call this the restricted adaptive variational quantum dynamics (RAVQDS).On the other hand, an adaptive process can still be executed at every time step provided an operator in the operator pool lowers the McLachlan's distance.In other words, no matter what time step we are at, an operator will be added if it can lower the distance (by a relative threshold).Such a protocol ensures that the McLachlan's distance will reach the lowest possible value possible and not be restricted at a fixed threshold.We therefor call this the unrestricted adaptive variational quantum dynamics or UAVQDS.
In the generic case, the RAVQDS algorithm is not applicable to the open system case due to the lack of a reasonable threshold, even when there are infinite operators in the pool.We present a sufficient condition under which the RAVQDS could be applied, assuming the operators in the pool are expressive enough, and defer the proof to Appendix G.A lower bound of the McLachlan distance square exists if Eq. (21) seems to be a very strong condition.However, it is trivially satisfied when solving a quantum trajectory equation Eq. ( 8) with Lindblad operators as Pauli operators.

Algorithms
In this section we present the pseudocode of our algorithms.Algorithm 1 and Algorithm 2 show the vectorization method and the trajectory method respectively.The adaptive procedure in both these algorithms can be either restricted or unrestricted, depending on whether Eqs.(21) are satisfied.The pseudocode for the unrestricted procedure in presented in Algorithm 3. In all the Algorithms, we use notations θ and A to denote the vector of the ansatz parameters (see Eq. ( 15)) and the corresponding ansatz operators.Additionally, Γ is the integrated norm shrinking factor given by Γ ≈ i dΓ i where where dΓ i is the norm shrinking factor of the ith step (Eq.( 19)).

Algorithm 1: Vectorization
Before presenting the simulation results, we would like to make a few comments on the algorithms used.In the vectorization method, we need to update θ and A at each time step and keep track of the norm e −Γ .The final vectorized state is e −Γ/2 |ϕ(t f )⟩.However, to measure any observable O in the unvectorized setting, we will need to evaluate where O † is the normalized and vectorized operator O † and ∥•∥ 1 denotes the trace norm.This Randomly pick a jump operator L i according to the probability mass function given by Eq. ( 12); quantity can be measured either directly or indirectly with the help of a synthesized unitary V which prepares O † , i.e.O † = V |0⟩ ⊗2N [32] (See Appendix H for details).
In the trajectory method, the total evolution is divided into intervals of a parameterized circuit followed by a jump operator.Each interval can be regarded as the state preparation circuit for the next one.As a consequence, the circuit depth of this approach increases with the number of jumps.Because we only need to measure the normalized state vector in this case, any measurement circuit described in Appendix E can be used.
Lastly, a threshold r which specifies the minimum amount reduction in the McLachlan distance needed to keep the adaptive procedure running is still necessary in the unrestricted algo- rithm.It can be either a additive threshold or a multiplicative threshold.In this paper, we will use the additive threshold as described in Algorithm 3, and refer it as the adaptive threshold.In addition, the algorithms (restricted or unrestricted ) require us to solve the linear equations given by Eq. ( 16), which can become illconditioned as the size of ansätze increases.To address this challenge, Tikhonov (L 2 ) regularization is applied when inverting M.Here λ is a small parameter to shift the diagonals of the M † M matrix.

Noiseless simulation
The QA simulation is carried out by evolving a system from t = 0 to t = t f under a time-dependent Hamiltonian of the form where A(t) and B(t) (known as the annealing schedule) are scalar functions which satisfy Here, we will only consider the linear schedule H D and H P are constant terms given by where N is the total number of qubits and Z i (X i ) is the Pauli Z (X) matrix on the ith qubit.
The couplings strength j i s are given by where n is the sector size.In this study, we fix the model parameters as w 1 = 1, w 2 = 0.5 and n = 1.
We focus on two types of Lindblad models, with the first one consisting of only Z i Lindblad operators applied to each qubit, given by: with rate γ i .It describes a continuous dephasing channel on each qubit.It is worth noting that the conditions outlined in Eqs.(21) are satisfied for the corresponding unravelled equation in this case.Therefore, the RAVQDS method can be applied for the trajectory method.We will refer to this model as the dephasing model going forward, and we will fix the values of dephasing rates as γ i = 0.01 in our simulations.The second model consists of two Lindblad operators on each qubit with different rate γ + i and γ − i .This model is often referred to as an empirical model for an incoherent energy transfer process and can also be rigorously derived from first principles.We will refer to this model as the amplitude damping model, and fix the values of the rates as γ + i = 0.04 and γ − i = 0.004 henceforth.It is worth noting that we choose the γ values based on two considerations: First, they should be small enough to ensure the weak coupling limit still holds; Second, they need to be large enough to make the open system effects non-negligible at the time scale of our simulation.An additional challenge arises when implementing the trajectory simulation of the amplitude damping model, because the jump operators L + i and L − i are non-unitary.Either midcircuit measurement or block encoding (unitary dilation) [54][55][56] based methods can be adopted to implement those non-unitary jumps (See Appendix I for details).
We conducted numerical experiments to compare different operator pools.The choice of our operator pool is a combination of the Hamiltonian pool [43] and the qubit-adapt pool [57].Initially, we utilize the smallest pool, namely the Hamiltonian pool, and incrementally enlarge it if the algorithm does not achieve the desired level of accuracy.While there can be multiple choices of operator pools, variational ansätze generated with qubit-adapt pools are much shallower [39,57].Starting from a pool consisting of all single qubit Pauli operators , we constructed three distinct operator pools by adding three different types of two-qubit Pauli operators to P s , defined as follows: where and N ′ denotes the number of physical qubits.We use the individual terms from Eqs. (30) as the operators in our pool.P 1 draws its inspiration from the Hamiltonian pool (the terms in the Hamiltonian appear in the pool only).We found the P 1 fails to recover the desired dynamics, which could be because a Hamiltonian pool is effective only for simulating unitary dynamics.Since the problem of open quantum system deals with non-unitary components as well, a Hamiltonian pool is insufficient in this case.While P 2 generates the desired results for the trajectory method, it fails in the vectorization method case, which could be attributed to the non-local term Lk ⊗ L k in the vectorized effective Hamiltonian (Eq.( 7)).Only P 3 proves effective in producing the desired results for both methods.We summarize the pool dependence of our method in Table 1.
When designing an operator pool, an important consideration is identifying the minimal com- plete pool (MCP), which minimizes the extra measurements needed for the adaptive process.MCP has been extensively studied in the case when the state vector is real [57,58], where an explicit MCP is suggested.However, for the problem we study in this paper (Eq.( 26)), the state vector is not limited to real vectors.As a result, the MCP proposed in the aforementioned references is are not directly applicable to our case.Although we can reformulate the problem such that the existing MCP results would apply, i.e., the state vector is always real, by expressing Eq. (1) in the Pauli basis, we have opted not to pursue this route in our current study.This is because such a representation would lead to an effective Hamiltonian comprising 4-body Pauli terms, potentially increasing circuit complexity (see Appendix C for a detailed discussion).Instead we provide an upper bound on the number of operators necessary for an MCP applicable to any state vector by selecting the operator pool such that it can generate a universal gate set.One such pool comprises every single Pauli operator and the ZX operator on adjacent qubits, leading to an MCP size upper bound of 4n − 1.
On the other hand, as indicated by [58] and by our numerical findings presented in Table 1 (P 2 is already a complete pool), having a complete pool does not guarantee optimal outcomes.Investigating strategies that enhance the algorithm's ability to converge towards accurate results represents a promising direction for further research.
In Figs. 2 and 3, we present the results from noiseless simulations of running UAVQDS with both the trajectory and vectorization methods on the dephasing and amplitude damping models.We compare the evolution of the energy ⟨H(t)⟩ = Tr{ρH(t)} and instantaneous eigenstate populations P ε i (t) = ⟨ε i (t)|H(t)|ε i (t)⟩ obtained from UAVQDS with the exact solution obtained from the Hamiltonian open quantum system toolkit (HOQST) package [59].Here, |ε i (t)⟩ represents the ith instantaneous eigenstate of the Hamiltonian, i.e., H(t) |ε i (t)⟩ = E i |ε i (t)⟩.The eigenstate populations are chosen as benchmarks because they are the quantities of interest in the original QA setting.The non-vanishing excitedstate populations indicate that the dynamics happens in the non-adiabatic regime, where quantum effects are expected to play a non-trivial role [60][61][62][63][64]. Since the primary focus of our work in not QA, we direct readers keen on further details to the references mentioned above.
While it is still unclear whether the eigenstate populations can be efficiently measured on a real quantum computer, they can be easily obtained in the simulation.We will focus only on the energy evolution when implementing the algorithm on the hardware.When running the simulations, we choose a stepsize of dt = 0.01 with a total evolution time t f = 10.Additionally, we set the Tikhonov regularization parameter λ in Eq. (23) to a fixed value of 10 −8 .Before proceeding, it's important to note that the instantaneous eigenstate populations depicted in Figs. 2 and 3 do not sum to one as they represent only the lowest 8 energy states.We have confirmed that the trace of the density matrix throughout the evolution remains unity.
The figures demonstrate that UAVQDS produces results in good agreement with the exact solver, as the dots representing the exact solutions overlap with the UAVQDS lines across all plots.It is noteworthy that compared to the trajectory method, the vectorization method requires double the number of qubits and a bigger operator pool size with non-nearest neighbor coupling.Hence, the ansatz there is more complex and more expensive to generate.For instance, r = 10 −4 was used to run the (N = 8) trajectory problem (with 8 physical qubits), whereas for the (N = 4) vectorized problem (also with 8 physical qubits), r had to be reduced to 10 −6 (10 −7 for the amplitude damping model) to ensure enough operators were included in the simulation (See Appendix K for more data).This implies that the effective Schrödinger equation generated by the vectorization method poses a greater challenge for the UAVQDS than the one produced by the trajectory method.This introduces an additional trade-off between the vectorization and trajectory methods, in addition to computational space (number of qubits) and execution time (number of trajectories).
We further investigate the error of the vectorization UAVQDS in Fig. 4. To quantify the error, we define an infidelity measure D as the tracenorm distance between the solution of UAVQDS and the exact solution at the end of the evolution, given by where ρ(t f ) and ρ t represent the solutions obtained from UAVQDS and the exact solver, respectively.As shown in Fig. 4a, by lowering the adaptive threshold, we increase the accuracy of the algorithm at the cost of more ansatz operators.However, this trend will not continue indefinitely because there are only a finite number of operators in the predefined pool.The relationship between the infidelity and total number of parameters in the ansätze is illustrated in Fig. 4b.The infidelity decreases approximately polynominally (linearly in log-log plot) as more operators are included in the ansätze until it hit a flat region.To further reduce the infidelity, a larger pool is required.Finally, we discuss the resource requirement for our algorithms.Two key factors to consider are the total number of operators and the number of multi-qubit operators in the ansatz, both of which are positively correlated with the final circuit depth required to execute the algorithm.To estimate the maximum number of CNOT gates needed for each time step, we start by calculating the number of CNOTs in the ansatz.The ansatz consists of unitaries that have the form e −iθP l , where P l is a Pauli word of length l.Thus, the total number of CNOT gates needed is given by l>1 2(N l − 1) where N l is the number of P l s in the ansatz.The scaling of the ansatz and estimated CNOT counts versus the system size N is shown in Figs.5a and 5b for both the trajectory and vectorization methods.Please refer to Appendix J for a demonstration on how the sizes of the ansätze grow as operators are adaptively added.The largest CNOT count in the figure is around 900 which is still far below the required CNOT count for any Trotterized algorithm using the same step size (with dt = 0.01, at least 1000 trotter step is needed).
In order to study scaling of parameters and resources, the operator counts are fitted to either a polynomial of the form aN b or an exponential of the form αe βN using least squares algorithm, and the results are summarized in Table.2. To compensate for the non-monotonic behavior of operator counts in the vectorization method, we add (0, 0) to the fitting data.Since the operator count profiles for different trajectories in the trajectory method vary, we compute fitting data based on the (a) mean, (b) median and (c) maximum value of the operator counts at each time step of the trajectories.Based on the coefficient of determination (R 2 ) for both the polynomial and exponential fittings, it is clear that a low degree polynomial better fits our data.However, more data points are needed if we want to rule out the exponential model with high confidence.The curves for polynomial fitting are displayed in Figs.5c and 5d.Based on the scaling data for the finite system sizes we studied, we do not observe signs of exponential scaling.AVQDS has previously shown polynomial scaling in terms of the number of parameters and CNOTs for other spin models [43].Confirming this observation for larger system in case of open quantum systems will be a topic for future studies.

Hardware results
To showcase our algorithm on a quantum computer that is currently available, we save the parameters of the time evolution that are computed classically using UAVQDS, and use these parameters to measure the energy at each time step by executing the circuits for Eq.(24) on a real quantum computer.We run our experiments on IBM's 27-qubit processor ibmq_kolkata with the  Falcon architecture.To ensure the algorithms run smoothly, it is necessary to carefully compile the circuit based on the specific quantum device be-ing used and to apply error suppression and mitigation strategies to obtain accurate and dependable results.
We note that we only simulate the dephasing model on the hardware because we want to avoid the additional complexity of mid-circuit measurement or reset operations required by simulating the amplitude damping model using the trajectory method (see Appendix I).Additionally, we compare the vectorization and trajectory method by selecting their respective N s such that they both use 4 physical qubits, meaning N = 2 for the trajectory method and N = 4 for the vectorization method.
Circuit generation-We utilize the Qiskit transpiler to generate multiple circuits with varying numbers of CNOTs due to the stochastic addition of swap gates, and select the circuit with the lowest number of CNOTs for our study.Direct measurement is used for the vectorization method, while indirect measurement via the Hadamard test is used for the trajectory method (see Appendix H).To synthesize the unitary required for measurements in the vectorization method, we use the Berkeley Quantum Synthesis Toolkit (BQSKit) [65,66].We do not perform additional recompilation besides the optimization supported by the Qiskit transpiler.
Error mitigation and post processing-We use the matrix-free measurement mitigation (M3) package to apply readout error mitigation in our study.This approach operates within a reduced subspace that is defined by the noisy input bitstrings requiring correction, making it scalable.The M3 package is described in references by [67] and is accessible through Qiskit.We also incorporate dynamical decoupling using Qiskit's standard tools, implementing periodic gate sequences to suppress cross-talk and system-environment coupling [68][69][70][71].Additionally, we apply an empirical error mitigation scheme based on the classical resolution enhancement (RE) technique, which has shown effectiveness in a different context [72], to the trajectory simulations.In order to avoid the ambiguity of the value of enhancement factor, the RE method uses a convergence criteria based on tracing over the measured probabilities of the non-ancillary qubits in a Hadamard test setting.Since the vectorization method uses direct measurement without any ancillary qubit, we were unable to apply RE in this case.It is important to note that all of the error mitigation strategies we applied are out-of-thebox solutions and are resource-efficient, making them suitable for real-world applications.
Results-We show the experimental results on the hardware in Fig. 6.The figure shows ⟨H⟩ versus evolution time t for vectorization (a) and trajectory (b) methods, respectively.It includes results obtained from ideal circuit simulation, quantum computer with and without error suppression or mitigation techniques, as well as the exact solver.The vectorization result is for a problem size of N = 2, i.e, 4 physical qubits.The trajectory method is for a problem size N = 4, which requires 4 physical qubits plus one additional ancilla qubit in order to measure ⟨H⟩ with the Hadamard test.
From Figs. 6a and 6b we observe good qualitative agreement between the exact solution and the results obtained from both the vectorization and trajectory methods.It is important to note that the error bars for the vectorization method and trajectory method represent distinct quantities.Specifically, the error bars for the vectorization method represent two times the standard deviation of 3 runs, while the error bars for the trajectory method represent two times the standard error of the mean over 17 trajectories.The performance of the vectorization method is particularly impressive because one would expect it to perform worse than the trajectory method, given that the N = 2 vectorization problem requires more operators in the ansatz (see Fig. 5) than the N = 4 trajectory problem and its ansätze include non-local 2-qubit Pauli operators.In addition, we expect the results to improve significantly with aggressive circuit recompilation since minimal fine-tuning was done in our hardware runs.The results from the trajectory method improve considerably after post-processing with the RE technique.However, it is important to note that we have only explored the RE method for a few limited cases and its generalizability as a potential error mitigation technique is still not fully understood.Summarizing, we see that the results from the IBM quantum computers are highly encouraging and suggest that simulating Lindblad equations is within the capability of NISQ devices.

Conclusion
We have introduced an adaptive variational quantum algorithm to simulate open quantum system dynamics described by a Lindblad equation.The method variationally solves either a vectorized or a stochastically unravelled Lindblad equation , with the ansatz built adaptively at every time step by minimizing McLachlan's distance in a greedy manner.We provide a benchmark of the algorithm's performance on both the ideal simulator and IBM's quantum processor for simulating the open-system quantum annealing dynamics, achieving good quantitative and qualitative agreement with the exact solution.Additionally, we perform a resource analysis on finite systems and find polynomial scaling, indicating the algorithm's potential to be extended to larger sys- tems.Based on our numerical findings, we conjecture that both algorithms will exhibit polynomial scaling in gate complexity, given the appropriate selection of an ansatz pool.Additionally, these algorithms demonstrate a space-time trade-off (number of trajectories versus number of qubits) that mirrors what we observed in their classical counterparts.When considering NISQera quantum devices, the trajectory method is more favorable because of two reasons.First, qubits are more valuable as a resource in the NISQ era.Secondly, the vectorization method presents an additional challenge due to its nonlocal effective Hamiltonian, which requires more complex ansätze.
Scaling our algorithm to sizes where meaningful simulations [4,14,15,73] can be conducted on NISQ hardware presents several challenges, with the foremost being the required circuit depth.Despite considerable advancements in enhancing the quality of quantum processors, decoherence still sets a practical limit on the achievable circuit depth.For instance, the most advanced experiments on IBM's 127-qubit chip have been constrained to circuit depths of less than 100 [74].In contrast, our algorithm's circuit depth surpasses this threshold even for prob-lems involving just a few qubits, as illustrated in Fig. 5, where the circuit depth closely follows the CNOT count.Besides waiting for better hardware, enhancements in both algorithmic design and software optimization are essential.For example, employing advanced recompilation techniques with heuristic algorithms can substantially reduce circuit depth [65]; implementing advanced error mitigation strategies can allow for longer circuit [22,75].Additionally, the training part of our algorithm, currently performed classically, will ultimately need to be transitioned to a quantum computer.This transition will benefit from identifying the MCP to reducing the measurement overhead of the training process and developing strategies that ensures convergence of the training [58] The authors also acknowledge the use of IBM Quantum services.The views expressed are those of the authors and do not reflect the official policy or position of IBM or the IBM Quantum team.
The supplementary code for this paper is available in the this GitHub repository.

A Frequently used symbols
For convenience, we provide a table summarizing frequently used symbols in the main text.
Coefficient matrix and vector of the linear equation generated by the variational principle Γ A classical register to record the norm of ψ(t) in the algorithm Total number of parameter in the ansätze at the end of evolution D Trace-norm distance between the solution of variation algorithm and the exact solution

B Derivation of the vectorized Lindblad equation
In this section we show how to derive the vectorized Lindblad equation from where the dissipative term The strategy is to apply the identity vec(ABC) = (C T ⊗ A) vec(B) to each term on the right-hand side of Eq. (32).The Hamiltonian part [H(t), ρ(t)] can be written in terms of The Lindblad part can be written in terms of Summing up every term in Eqs.(34) and (35), the effective Hamiltonian can be obtained as C Vectorization in the Pauli basis An alternative approach to vectorizing the ME involves expressing the density matrix in the basis of Pauli matrices and representing the superoperators through the Pauli transfer matrices (PTMs).However, this method leads to a notable increase in circuit complexity for two primary reasons.
• Vectorization can transform a single-qubit pure state into a more complex state, such as a 2-qubit Bell state.For instance, consider the state |0⟩⟨0|.In the Pauli basis arranged in lexicographic order, this state, when vectorized, would be represented as (|00⟩ + |11⟩)/ √ 2 (up to a constant), necessitating an additional layer of 2-qubit gates for state preparation.
• When employing the column-stacking method for vectorization, if every term in the Hamiltonian is 2-local and each Lindblad operator is 1-local, then each term in the resultant effective Hamiltonian (as per Eq. ( 36)) remains at most 2-local.For example, the PTM for the superoperator [ZZ, •] is To implement the Hamiltonian pool [43] would require an arbitrary angle 4-qubit Pauli rotation, substantially increasing the number of required CNOT gates.
Given these considerations, particularly the worse error rates associated with 2-qubit gates on contemporary quantum devices, we opted against this vectorization approach to reduce the circuit depth.The potential for a more efficient implementation of this approach remains an intriguing area for future research.

D Derivation of the evolution equation
In this appendix, we derive the evolution equation of the variational parameters for the effective Schrödinger equation where the dot symbol denotes the time derivative.Before proceeding, it is worth mentioning that the effective Hamiltonian in Eq. (37) is not necessarily Hermitian.The ansatz state is parameterized by where A µ and |ψ R ⟩ are the ansatz operators and the reference state, respectively.It is more convenient to derive the evolution equation of θ(t) using the density matrix representation to avoid the problem of a time dependent global phase [50].Let W (θ) ≡ |ϕ(θ)⟩⟨ϕ(θ)|, the equation of motion is given by where we make use of Eq. (37) in going from Eq. (39a) to Eq. (39b).Applying the McLachlan's variational principle to Eq. (39), we have where Based on Eq. (39a), Ẇ is Hermitian, i.e., Ẇ = Ẇ † .The first term Tr Ẇ † Ẇ on the RHS of Eq. ( 41) can be simplified by the chain rule where M is a matrix with elements Because θ µ are real parameters, M is a real symmetric matrix.
Similarly, the second term on the RHS of Eq. (41) can be simplified as where V is a vector with elements and ⟨H eff ⟩ = ⟨ϕ|H eff |ϕ⟩.
Before calculating the last term on the RHS of Eq. (41), we split the effective Hamiltonian into a Hermitian and an anti-Hermitian parts Then the last term becomes The variational principle (Eq.( 40)) then yields Finally, the evolution equation of the variational parameters is given by a linear equation E Measurement circuit for M and V In this appendix, we will discuss the quantum circuit used to measure the parameter in the variational equation of motion M θ = V, where and assuming that each ansatz operator A µ is a Pauli operator.We only present the resulting circuits here and encourage interested readers to refer to [53] for further details.First, the term ⟨H eff ⟩ can be evaluated by measuring all the Pauli strings P i that make up the Hamiltonian.The quantum circuit for both direct and indirect measurement (Hadamard test) are shown in Fig. E.1, where we use notation U ⃗ θ to denote the the unitary generated by the variational circuit and H S † to denote the optional Hadamard (Hadamard-phase) gate to rotate the basis of P i to z direction.Second, the term ⟨ϕ| ∂|ϕ⟩ ∂θµ and ⟨ϕ| H eff ∂|ϕ⟩ ∂θµ can be evaluated by using a generalized Hadamard test or the corresponding direct measurement circuit.The main observation is that, using the term ⟨ϕ| P i ∂|ϕ⟩ ∂θµ can be simplified as where U µ:1 = µ j=1 e −iθ j A j and U k:µ+1 = k l=µ+1 e −iθ l A l .Eq. ( 53) can be evaluated using the generalized Hadamard test (Fig. E.2a) or the corresponding direct measure circuit (Figs.E.2b and E.2c).Finally, the metric tensor, given by ∂⟨ϕ| ∂θµ ∂|ϕ⟩ ∂θν , can also be evaluated using circuits similar to those discussed above.We note that this expression can be simplified as

F Evolution of the state-vector norm
In this section, we show how to track the evolution of the state-vector norm by measuring the anti-Hermitian component of the effective Hamiltonian at every time step.We start with the Schrödinger equation with ψ(t) being the unnormalized state vector and |ψ(t)⟩ being the normalized state vector.The effective Hamiltonian has both the Hermitian H e and anti-Hermitian parts −iH a , i.e., H eff = H e −iH a .We also use the notation |ϕ(t)⟩ for the state generated by our variational circuit.The evolution of ψ(t) ψ(t) can be calculated as The solution to the above equation is In the above equation, e −2 t 0 ⟨ψ(τ )|Ha|ψ(τ )⟩dτ is monotonically decreasing.To see this, we first examine the case of the unravelled Lindblad equation, where k L k is positive semi-definite since each term in the summation is quadratic.As a result, t 0 ⟨ψ(τ )|H a |ψ(τ )⟩ dτ is always a positive value.We then turn to the case of the vectorized Lindblad equation.Because the evolution operator of the original Lindblad equation is a contraction map, i.e., and the linear isometry preserves the trace norm and L 2 norm, we have For a small time interval dt thus we can keep track of the norm by measuring H a at each time step.

G Lower bound of the McLachlan's distance
In this Appendix, we derive a lower bound of the McLachlan's distance (Eq.( 41)).First, we note from Eq. (38) that the ansatz state is parameterized by a unitary A reasonable assumption we make is that dU (θ(t))/dt exists (well-defined and finite).With this assumption, we can define an effective Hamiltonian parameterized by the derivative of the ansatz parameter where we omit the t-dependence of each quantity in the above equation.It is straightforward to check the following relation is satisfied Then the derivative of the ansatz state is subject to the following equation of motion and the McLachlan's distance becomes where we denote ∆H ≡ H − H e in going from Eq. (65b) to Eq. (65c).When [H a , ∆H] = 0, a lower bound of D is given by where the equality is achieved when ∆H = 0, which means the optimal θ makes H θ in Eq. (63) equal the Hermitian part of the effective Hamiltonian H where Eq. (67b) guarantees H θ , H a = 0 .

H Measure observable O in the vectorization method
In this section we describe the circuit used to measure the observable O for the vectorization method based on Ref. [32].The target is to evaluate the quantity O † ϕ(t) where O † is the vectorized operator O † from a smaller Hilbert space.The main observation is that the parameterized state |ϕ(t)⟩ can be generated from |0⟩ (all zero state) with one unitary where U R prepares the reference state.
where x 1 and y 1 are bit strings for 0 to 2 N − 1 and O x 1 y 1 is the x 1 , y 1 element of the operator.V can be synthesized using quantum state preparation algorithm [66,76].Because O is usually a local observable, Eq. ( 69) can be further simplified to avoid any exponential cost [32].In the examples considered in the main text, the reference state is chosen as |ψ R ⟩ = |+⟩ ⊗2N , as a result, U R can be implemented using Hadamard gates H ⊗2N .
In our experiments, we used the Berkeley Quantum Synthesis Toolkit (BQSKit) [65,66] to synthesize the unitary V and used the direct measurement protocol shown in Fig. H.1 to evaluate O † ϕ(t) .The Hadamard test is not currently practical on hardware because it requires implementing V † U ⃗ θ H ⊗2N as a control unitary.We can obtain the final result by taking the square root of the all-0 string frequency since O † ϕ(t) = ⟨0|V † U ⃗ θ H ⊗2N |0⟩.When taking the square root, a sign needs to be assigned; however, we can always shift an observable by the identity matrix such that all its eigenvalues have the same sign.

I Implementing L − and L + jump operators
In this appendix, we describe two methods for applying the L + and L − jump operators on a quantum computer.Without loss of generality, we present only the derivation for L + ≡ (X + iY )/2, and the procedure is easily extended to L − ≡ (X − iY )/2.The key observation is that the L + jump only occurs when its probability p + ∝ ⟨ϕ|L − L + |ϕ⟩ is greater than zero, meaning |ϕ⟩ ̸ = |0⟩.Under this condition, the action of the jump operator L + |ϕ⟩ / ⟨ϕ|L − L + |ϕ⟩ is equivalent to a quantum eraser channel E(ρ) = |0⟩⟨0|.There are two ways of implementing such a channel on a quantum computer: mid-circuit measurement or block encoding (unitary dilation) [54][55][56].The former requires us to perform a mid-circuit measurement on the qubit and apply a X gate conditioned on the result 1.The latter requires an additional ancilla qubit prepared in |0⟩ state.Then we apply a CNOT gate with the ancilla as the target qubit, measure the ancilla and post select the runs with measurement result 0. The circuits for both of these methods are shown in

J Growth of ansatz size
In this appendix, we demonstrate how the sizes of ansätze grow as we adaptively add operators to them.The results are obtained from the examples with N = 8 showcased in Fig. 5a and N = 4 in Fig. 5b in the main text.In Fig. J.1a, we present examples of three different trajectories where 0, 1, and 2 jump events occurred during the evolution.As we can see, the adaptive procedure quickly refills the ansätze after they are reset by the jumps.Consequently, the accumulated ansatz size becomes larger when more jumps occur.In

K The vectorization method with different operator pools and error threshold
In this appendix, we display the average energy ⟨H(t)⟩ versus evolution time t for both the dephasing and amplitude damping models obtained using the vectorized UAVQDS with different choices of operator pools and error thresholds.The purpose is to provide visualization of how the error behaves with respect to various operator pools and error thresholds.

L Resolution enhancement technique
Our method for resolution enhancement relies on the assumption that the bitstring data produced during a noisy experiment somewhat conforms to the expected probability distribution of bitstrings.
In simpler terms, even though noise and measurement errors cause the histogram of bitstrings to become less distinct, the underlying pattern of the histogram remains intact.To align the bitstring distribution with the ideal one, we apply resolution enhancement techniques commonly employed in image processing [77].The corrected bistring distribution data is then used to compute the expectation values.While the details of the method have been already described in [72], for the sake of completeness we briefly discuss the method here as well as provide the specific steps we have undertaken for our current work.
In our work, we measure a single Pauli term using Hadamard test circuit and later add multiple expectation values to find the energy expectation value ⟨H(t)⟩.Hadamard test requires one ancilla qubit and the expectation value can be obtained by just measuring the ancilla qubit.We, however, measure all physical qubits (including the ancilla) and then trace over the rest to obtain the state of the ancilla.Such additional measurements also help getting rid of certain unphysical states that arise in quantum experiments, e.g, keeping track of fermion number conservation [78].In our work such screening were not necessary.Writing explicitly, if the measured wavefunction of the full system is given by |ψ⟩ = j α j |b n ..b 1 ⟩ |0⟩ a + β j |b n ..b 1 ⟩ |1⟩ a , where b j s are the measured state of of the qubits, after applying our RE method, we will obtain a corrected wavefunction |ψ c ⟩ = j α (c) Calling y j = |α j | 2 as the frequency of the noisy data of j-th bitstring, resolution enhanced frequency is obtained using r j = y j − k 2 y ′′ j , were r j = ( α 2 ) is the reformed frequency and y ′′ j is the second derivative of the noisy data w.r.t the decimal representation of the bitstrings.The parameter k 2 can be modified to tune the resolution of the final data.The weighting factor k 2 can be chosen based on what gives the best trade-off between resolution enhancement, signal-to-noise degradation, and baseline flatness.The optimum choice depends upon the width, shape, and digitization interval of the the bitstrings.After obtaining r j we identify the bitstrings nearest to the peaks from the resolution enhanced data.We then switch back to the binary representation of the bitstrings and replace y j s by the r j s.
In order to avoid the ambiguity of the optimal value of k 2 , we iterate over several values of k 2 and calculate the probability p 0 of the ancilla qubit for each of them.We continue iterating until p 0

Figure 1 :
Figure 1: An illustration of how a non-unitary evolution is simulated using an unitary.Assume at time t the state is a pure state |ϕ(t)⟩ on the Bloch sphere (blue arrow).After an infinitesimal step dt, an evolution subject to a non-Hermitian effective Hamiltonian will bring the state inside the Bloch sphere (red arrow) by shrinking its norm ψ ψ < 1.The goal of the variational algorithm is to evolve |ϕ(t)⟩ to the closest state to ψ(t + dt) on the Bloch sphere (red cross), i.e., |ϕ(t + dt)⟩ = |ψ(t + dt)⟩, while recording ψ(t + dt) as classical information.

Figure 2 :
Figure 2: Noiseless simulation using the trajectory method and UAVQDS for N = 8 (8 physical qubits).Energy and eigenstate population evolution using the dephasing model ((a) and (b)) and the amplitude damping model ((c) and (d)) are shown respectively.The UAVQDS results are obtained from 1000 trajectories and the ribbons represent the 2σ (σ is the standard error of the mean) error bar for the trajectory average.The adaptive threshold was set to r = 10 −4 and the time step size was set to dt = 0.01.We used the operator pool P 2 .

Figure 3 :
Figure 3: Noiseless simulation using the vectorization method and UAVQDS for N = 4 (8 physical qubit).Energy and eigenstate population evolution using the dephasing model ((a) and (b)) and the amplitude damping model ((c) and (d)) are shown respectively.For the dephasing model, the adaptive threshold was set to r = 10 −6 , and for the amplitude damping model, it was set to r = 10 −7 .The time step size was fixed at dt = 0.01.The operator pool P 3 was used for these simulations.

Figure 4 :
Figure 4: Infidelity of UAVQDS versus adaptive threshold and total number of parameters in the ansatz.Panel (a) shows the total number of parameter in the ansätze at the end of evolution (blue dots) and the infidelity (red crosses) versus the adaptive threshold r.The operator pool employed includes all 2-qubit Paulis (P 3 ).Panel (b) shows the infidelity versus the total number of parameters in the ansätze for pools with either neighboring (blue dots) or all (orange dots) 2-qubit Paulis (P 2 and P 3 , respectively).The dashed line in (b) represents a fitted line using all data from the orange dots except the last one.The results were obtained using the setup from Fig.3a.

Figure 5 :
Figure 5: Resource requirements for different system sizes.(a) Ansatz and estimated CNOT counts for the trajectory method.(b) Ansatz and estimated CNOT counts for the vectorization method.(c) Ansatz counts fitted to a model of aN b for the trajectory method, with box plots of 1000 trajectories.(d) Ansatz counts fitted to a model of aN b for the vectorization method.In all the panels, N p represents the ansatz count, while CNOTs indicates the estimated CNOT count.The box plots in panels (a) and (c) are obtained from 1000 trajectories, and panels (b) and (d) display resource counts for two adaptive threshold choices.The results were obtained using the configuration described in Fig. 3a.

Figure 6 :
Figure 6: Energy expectation value obtained by executing the parameterized circuits on ibmq_kolkata.Panels (a) and (b) show the results for the vectorization method (N = 2) and trajectory method (N = 4), respectively.The noiseless, noisy, and exact curves represent the results from ideal circuit simulation, quantum computer without any error suppression or mitigation techniques and exact solver, respectively.Other curves show the results obtained from the quantum computer with corresponding error suppression or mitigation techniques: measurement error mitigation (M3), dynamical decoupling (DD) and resolution enhancement (RE).The error bars in panel (a) represent two times the standard deviation from 3 runs, and the error bars in panel (b) represent two times the standard error of the mean from 17 trajectories.All curves from IBM quantum computers are obtained with 100000 shots.Note that we do not have a noisy curve for the trajectory method because running trajectories on the real quantum computer is resource-intensive.

Table A. 1 :
List of symbols and their definitions used in the main text.

Figure E. 1 :
Figure E.1: Quantum circuit for evaluating ⟨ϕ|P i |ϕ⟩ using (a) Hadamard test; (b) direct measurements.U ⃗ θ denotes the unitary generated by the variational circuit k µ=1 e −iθµ(t)Aµ .H S † denotes the optional Hadamard (Hadmard-phase) gate to rotate the basis of P i to z direction.The gates within the blue block implement a projective measurement onto the eigenbasis of P i .

Figure E. 2 :
Figure E.2: Quantum circuits for evaluating ⟨ϕ| P i ∂|ϕ⟩ ∂θµ using (a) generalized Hadamard test; (b) and (c) direct measurements.The generalized Hadamard test circuit uses a binary integer b ∈ {0, 1} as an input.If b = 0, the circuit measures the real part of ⟨ψ R |U † µ:1 U † k:µ+1 P i U k:µ+1 A µ U µ:1 |ψ R ⟩, and if b = 1, it measures the imaginary part.For direct measurement, circuit (b) measures the real part of the targeted quantity, while circuit (c) measures the imaginary part.

Figure H. 1 :
Figure H.1: Direct measurement scheme to evaluate O † ϕ(t) .The final result can be obtained by taking the square root of the all-0 string frequency.

Figure I. 1 :
Figure I.1:Quantum circuit for implementing an eraser channel E(ρ) = |0⟩⟨0| using (a) mid-circuit measurement; (b) block encoding.The measurement symbol in (a) means a conditional gate based on the measurement result.We only apply X gate when the result is 1.The arrow pointing to |0⟩ in (b) represents the post-selection.We only keep the runs where 0 is returned by the measurement.

Figure J. 1 :
Figure J.1: Number of ansatz parameters required for different evolution times in both the trajectory and vectorization methods.In (a), the numbers of parameters are plotted for three different trajectories, with vertical lines indicating jump events where the ansätze are reset.In (b), the numbers parameters for the vectorization method is shown for different error thresholds, with the turquoise dashed line indicating the total number of distinct operators in the operator pool.

Figure K. 1 :
Figure K.1:The energy evolution of the vectorization UAVQDS for different Lindblad models, error thresholds and operator pools.Subplots (a) and (b) show the results obtained using a all-to-all 2-qubit Pauli pool (P 3 in the main text) for the dephasing and amplitude damping models, respectively.Subplots (c) and (d) show the results obtained using a neighboring 2-qubit Pauli pool (P 2 in the main text) for the dephasing and amplitude damping models, respectively.

Figure L. 1 : 2 )
Figure L.1: Effects of resolution enhancement on bitstring distribution.Histogram of bitstrings of noisy simulation (a) before and (b) after applying resolution enhancement

Table 2 :
A summary of fitting results.The ansatz counts shown in Fig.5are fitted to either a polynomial of the form aN b or an exponential of the form αe βN using a least squares algorithm.The table reports the coefficient of determination (R 2 ) for both fittings and the exponent of the resulting polynomial.
Hermitian part of the effective Hamiltonian H eff H a −iH a is the anti-Hermitian part of the effective Hamiltonian H eff D Square of the McLachlan distance used in the variational principle ψ(t) Unnormalized state vector solution for the effective Schrödinger equation with H eff |ψ(t)⟩ Normalized state vector solution for the effective Schrödinger equation H eff If we can construct another unitary that prepares O † , i.e.O † = V |0⟩, then O † ϕ(t) = ⟨0|V † U (θ)U R |0⟩ can be measured either directly or indirectly (Hadamard test).Assume the operator O has dimension 2 N × 2 N , an explicit form of O † is given by