Adaptive variational simulation for open quantum systems

Huo Chen
Computational Research Division, Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA

Niladri Gomes

H.C. and N.G. contributed equally.
Computational Research Division, Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA

Siyuan Niu

Wibe Albert de Jong
Computational Research Division, Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA

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 ansätze 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 <cit.>. Recently there has been a surge of interest in developing efficient simulation algorithms for open quantum systems <cit.>. On one hand, the emergence of quantum technologies for building artificial quantum systems for computing <cit.>, sensing <cit.>, light harvesting <cit.>, 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 <cit.>, relaxation dynamics of many-body quantum system <cit.>, charge and energy transfer dynamics in molecular systems <cit.>, 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 cloud-accessible quantum computers <cit.> for open-system 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 <cit.>.
Based on the underlying techniques, existing quantum algorithms for simulating the Lindblad dynamics fall into five categories: unitary dilation <cit.>, variational simulation <cit.>, quantum imaginary-time evolution (QITE) <cit.>, Monte Carlo <cit.> and qubit-bath engineering <cit.>. Among them the variational simulation and QITE require less circuit depth and are considered more NISQ (noisy intermediate-scale quantum) <cit.> friendly.

Variational algorithms have been shown to be capable for solving real and imaginary time evolutions, non-Hermitian quantum mechanics for non-equilibrium systems, open quantum systems and  general first-order differential equations <cit.>. In this paper, we present a compact approach for solving the Lindblad equation using a time-dependent adaptive variational method. Our strategy to generate NISQ-friendly ansätze is built upon adaptive variational quantum dynamics simulations (AVQDS) <cit.>. 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. To facilitate this process, we have developed an unrestricted adaptive protocol that addresses the challenge of establishing a reasonable 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 <cit.>.
QA is a quantum meta-heuristic algorithm that is commonly implemented on analog quantum hardware. Although our algorithm is specifically designed for gate-based quantum computers, simulating QA hardware highlights its applicability to a wider range of use cases, such as analyzing artificial quantum systems.

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.

The structure of this paper is as follows. In Sec. <ref>, we introduce the Lindblad equation and describe our quantum algorithm,  unrestricted adaptive variational
quantum dynamics (UAVQDS), for solving it. In Sec. <ref> we report the performance and resource requirement of our algorithm based on the results from the noiseless simulator.
In Sec. <ref>, we report the performance of our algorithm on IBM quantum computers. We conclude in Sec. <ref>, and present additional technical details and derivations in the Appendices.

§ METHOD

§.§ 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

\begin{equation}
\label{eq:Lindblad_eq}
\frac{d}{dt}\rho\pqty{t} = -i\comm{H\pqty{t}}{\rho\pqty{t}}+\Lag\bqty{\rho\pqty{t}} \ ,
\end{equation}

where $\rho\pqty{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 $\Lag\bqty{\rho\pqty{t}}$ which is given by the form,

\begin{equation}
\Lag\bqty{\rho\pqty{t}} = \sum_{k=1}^{K} \gamma_k \pqty{L_k \rho\pqty{t} L^\dagger_k -\frac{1}{2}\acomm{L^\dagger_k L_k}{\rho\pqty{t}}} \ .
\end{equation}

The operators $L_k$ are known as the Lindblad operators, and are taken to be dimensionless. Since the Lindblad equation is a dissipative model, the rate of dissipation or decay of the system over time is modeled by varying $\gamma_k$. The $\gamma_k$ are non-negative quantities with dimension of inverse time. Usually it is convenient to absorb $\gamma_k$ into the definition of $L_k$,

\begin{equation}
\sqrt{\gamma_k} L_k \mapsto L_k \ .
\end{equation}

Throughout our discussions, we have adopted the convention of setting $\hbar = 1$.

Because a quantum computer only evolves a pure state under unitaries, the density matrix equation (Eq. (<ref>)) cannot be directly solved by a gate-based quantum computer.   This challenge can be addressed by two different approaches: 1) vectorizing <cit.> or 2) stochastic unravelling of the Lindblad ME <cit.>. We briefly introduce 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 $\rho$ on top of one another:

\begin{equation}
\mathrm{vec}\pqty{\rho} \equiv \bqty{\rho_{11}, \rho_{12}\cdots\rho_{21},\rho_{22}\cdots\rho_{\mathrm{D}\mathrm{D}}}^\mathrm{T} \ ,
\end{equation}

where $\mathrm{D}\times \mathrm{D}$ is the dimension of $\rho$. We will also use the symbol $\rvert \rho \rrangle$ interchangeably with $\mathrm{vec}\pqty{\rho}$. The map $\mathrm{vec}$ is a linear isometry between the $\mathrm{D}\times \mathrm{D}$ Hilbert–Schmidt space of $\rho$ and the $\mathrm{D}^2$ Hilbert space of $\rvert \rho \rrangle$, with the preserved norms being the trace norm and $L_2$ norm

\begin{equation}
\label{eq:iso_trace}
\sqrt{\Tr\pqty{\rho^\dagger\rho}} = \sqrt{\llangle \rho \vert\rho \rrangle} \ .
\end{equation}

Using the identity $\operatorname{vec}(A B C)=\left(C^{\mathrm{T}} \otimes A\right) \operatorname{vec}(B)$, Eq. (<ref>) can be rewritten as

\begin{equation}
\frac{d}{dt}\kett{\rho\pqty{t}} = -iH_{\mathrm{eff}}\pqty{t} \kett{\rho\pqty{t}} \ ,
\end{equation}

where $H_{\mathrm{eff}}$ is an effective non-Hermitian Hamiltonian given by (see Appendix <ref> for the derivation)

\begin{align}
\label{eq:Heff_vec}
H_\mathrm{eff} &= I\otimes H - H^T\otimes I \notag\\
&+i\sum_k \bqty{L_k^*\otimes L_k-\frac{1}{2}\pqty{I\otimes L_k^{\dagger} L_k + L_k^{T} L^*_k\otimes I}} \ .
\end{align}

We will refer to this approach as the vectorization method for the remainder of the paper.

Stochastic unravelling –
The Lindblad master equation can be unravelled using quantum trajectory method <cit.>. 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 (different from Eq. (<ref>))

\begin{equation}
\label{eq:unravelled_equation}
\frac{d}{dt}\ket{\tilde{\psi}\pqty{t}} = -i H_\mathrm{eff}\ket{\tilde{\psi}\pqty{t}} \ ,
\end{equation}

where $\ket{\tilde{\psi}\pqty{t}}$ is the unnormalized state vector and the effective Hamiltonian is given by

\begin{equation}
H_\mathrm{eff} = H\pqty{t} -\frac{i}{2}\sum_{k=1}^{K} L^\dagger_k L_k \ .
\end{equation}

For an infinitesimal time step from $t$ to $t+dt$, there are two possible evolutions for $\ket{\tilde{\psi}\pqty{t}}$: either the deterministic evolution subject to Eq. (<ref>) (with probability $1-dp$) or a jump occurring (with probability $dp$). The probability of a jump is given by

\begin{equation}
dp = \sum_{k=1}^K \expval{L_k^\dagger L_k} dt \ .
\end{equation}

Furthermore, if a jump happens, the unnormalized state is updated as

\begin{equation}
\ket{\tilde{\psi}\pqty{t+dt}} = L_i\ket{\tilde{\psi}\pqty{t}}/\sqrt{\mel{\tilde{\psi}\pqty{t}}{L^\dagger_i L_i}{\tilde{\psi}\pqty{t}}}
\end{equation}

where $L_i$ is randomly picked from $\Bqty{L_k}_{k=1}^K$ with probability

\begin{equation}
\label{eq:jump_p}
p_i =     \mel{\tilde{\psi}\pqty{t}}{L^\dagger_i L_i}{\tilde{\psi}\pqty{t}}/\sum_{k=1}^K\mel{\tilde{\psi}\pqty{t}}{L^\dagger_k L_k}{\tilde{\psi}\pqty{t}} \ .
\end{equation}

Let us denote the normalized state vector by $\ket{\psi_j\pqty{t}}=\ket{\tilde{\psi_j}\pqty{t}} / \sqrt{\braket{\tilde{\psi_j}\pqty{t}}{\tilde{\psi_j}\pqty{t}}}$ where the index $j$ represent the $j$th trajectory. The density matrix solution to Eq. (<ref>) can be constructed by
$\rho\pqty{t}=\frac{1}{n}\sum_{j=1}^n \dyad{\psi_j\pqty{t}}{\psi_j\pqty{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 <cit.>.
Without loss of generality, we split the effective Hamiltonian into a Hermitian and an anti-Hermitian parts and denotes them by $H_\mathrm{e}$ and $H_\mathrm{a}$ throughout this paper

H_eff = H_eff+H^†_eff/2 -i iH_eff-H^†_eff/2 ≡H_e - i H_a  .
Inspired by AVQDS <cit.>, 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 parameterised circuits

\begin{equation}
\label{eq:ansatz}
\ket{\phi\pqty{t}} = \prod_{\mu=1}^{k} e^{-i\theta_\mu\pqty{t} A_\mu} \ket{\psi_\mathrm{R}} \ ,
\end{equation}

where $e^{-i\theta_\mu\pqty{t} A_\mu}$ is the $\mu$th layer of circuit controlled by the real
parameter $\theta_\mu\pqty{t}$, and $\ket{\psi_\mathrm{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., $\theta_\mu\pqty{t}$ in Eq. (<ref>). For adaptive protocols, the ansatz operators $A_\mu$ are adaptively added from an ansatz pool at runtime.

Our algorithm is based on McLachlan variational
principle <cit.>

δd|ϕt⟩/dt+iH_eff|ϕt⟩^2 =0  ,
where $\bth\pqty{t}$ denotes $\bqty{\theta_1\pqty{t}, \theta_2\pqty{t}, \cdots \theta_k\pqty{t}}$. The norm in Eq. (<ref>) is known as the McLachlan distance and we denote its square by $\mathcal{D}$. The corresponding evolution equation of the variational parameters is given by

=  ,
where $\bM$ is a matrix with elements

_μν =2Re⟨ϕ|θ_μ|ϕ⟩θ_ν +⟨ϕ||ϕ⟩θ_μ ⟨ϕ||ϕ⟩θ_ν  ,
$\bV$ is a vector with elements

_μ= 2Im⟨ϕ||ϕ⟩θ_μ+⟨ϕ|θ_μ|ϕ⟩  ,
where $\expval{\Heff} = \mel{\phi}{\Heff}{\phi}$. Detailed derivation of the evolution equation can be found in the Appendix <ref>. Assuming every ansatz operator is a single Pauli string, each term in Eqs. (<ref>) and (<ref>) involving the derivatives can be evaluated using a quantum computer by either direct or indirect measurements <cit.> (See Appendix <ref> 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. (<ref>)) 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 <ref>, for an infinitesimal time step $dt$, the norm of $\ket{\tilde{\psi}\pqty{t}}$ shrinks according to

\begin{equation}
\label{eq:norm_evo}
\braket{\tilde\psi\pqty{t+dt}} \approx e^{- d \Gamma } \braket{\tilde\psi\pqty{t}} \ ,
\end{equation}

where $d \Gamma = 2\mel{\psi\pqty{t}}{H_\mathrm{a}}{\psi\pqty{t}} dt$, which is equivalent to $\mel{\phi}{H_\mathrm{a}}{\phi}$, can be measured using the same circuit that is used to measure the effective Hamiltonian. In Fig. <ref>, we provided a single qubit example using Bloch sphere to demonstrate of how our algorithm works.

An illustration of how a non-unitary evolution is simulated using an unitary. Assume at time $t$ the state is a pure state $\ket{\phi\pqty{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 $\braket{\tilde{\psi}}{\tilde{\psi}} < 1$. The goal of the variational algorithm is to evolve $\ket{\phi\pqty{t}}$ to the closest state to $\ket{\tilde{\psi}\pqty{t+dt}}$ on the Bloch sphere (red cross),
i.e., $\ket{\phi\pqty{t+dt}}= \ket{\psi\pqty{t+dt}} $,
while recording $\norm{\ket{\tilde{\psi}\pqty{t+dt}}}$ as classical information.

§.§ Restricted and unrestricted adaptive protocol

Unlike existing variational methods <cit.> 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 <cit.> 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 presented by in an adaptive protocol is that an attainable lower bound for the McLachlan distance is not known 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 <ref>.
A lower bound of the McLachlan distance square exists
𝒟 ≥2^2 + 2^2

\begin{align}
\comm{\He}{\Ha} &= 0 \\
\comm{A_\mu}{\Ha} &= 0 \ , \forall \mu \in \bqty{1, k} \ .
\end{align}

Eq. (<ref>) seems to be a very strong condition. However, it is trivially satisfied when solving a quantum trajectory equation Eq. (<ref>) with Lindblad operators as Pauli operators.

§.§ Algorithms

In this section we present the pseudocode of our algorithms. Algorithm <ref> and Algorithm <ref> 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. (<ref>) are satisfied. The pseudocode for the unrestricted procedure in presented in Algorithm <ref>. In all the Algorithms, we use notations $\bth$ and $\mathbf{A}$ to denote the vector of the ansatz parameters (see Eq. (<ref>)) and the corresponding ansatz operators. Additionally, $\Gamma$ is the integrated norm shrinking factor given by $\Gamma \approx \sum_i d\Gamma_i$ where where $d\Gamma_i$ is the norm shrinking factor of the $i$th step (Eq. (<ref>)).

$ dt $*[r]Step size
$\ket{\psi(0)}$ *[r]Initial state
$\bth,\ \mathbf{A},\ \Gamma$
$t \gets 0$
$\Gamma \gets 0 $ *[r]$\Gamma$ record the norm of $\ket{\tilde\psi\pqty{t}}$
$\bth, \mathbf{A} \gets []^T$*[r]Start with empty ansatz
$t < t_f$

update $\bth,\ \mathbf{A}, \ \dot{\bth}$ with the adaptive procedure
$\bth \gets \bth + \dot{\bth} dt$
$\Gamma \gets \Gamma + 2\expval{H_\mathrm{a}} dt$
$t \gets t+dt$

$dt$ *[r]Step size
$\ket{\psi(0)}$ *[r]Initial state
$\bth,\ \mathbf{A},\ \ket{\psi_\mathrm{R}}$

$t \gets 0$
$\Gamma \gets 0 $ *[r]$\Gamma$ record the norm of $\ket{\tilde\psi\pqty{t}}$
$\bth, \mathbf{A} \gets []^T$*[r]Start with empty ansatz
Generate a random number $q\in [0, 1)$
$t < t_f$
$e^{-\Gamma}\ge q$

update $\bth,\ \mathbf{A}, \ \dot{\bth}$ with the adaptive procedure
$\bth \gets \bth + \dot{\bth} dt$
$\Gamma \gets \Gamma + 2\expval{H_\mathrm{a}} dt$

Randomly pick a jump operator $L_i$ according to the probability mass function given by Eq. (<ref>)
$\ket{\psi_\mathrm{R}}\gets L_i\ket{\phi}/\norm{L_i\ket{\phi}}$
$\bth, \mathbf{A} \gets []^T$*[r]Reset the ansatz
$\Gamma \gets 0$
Generate a random number $q\in [0, 1)$

$t \gets t+dt$

Unrestricted adaptive procedure
$\ket{\phi}$ *[r]Parameterized state
$r$ *[r]Relative threshold
$\bth,\ \mathbf{A},\ \dot{\bth}$

Measure $\bM$ and $\bV$ in Eq. <ref>
$\dot\bth \gets \bM \backslash \bV$ *[r]We use Tikhonov regularization for the pseudoinverse
$\mathcal{D}'\gets\dot{\bth}^T \bM \dot{\bth}-2 \bV^T \dot{\bth}$
$\mathcal{D}' \ge \mathcal{D}-r$
$\mathcal{D} \gets \mathcal{D}'$
$\bM' \gets \bM, \ \bV' \gets \bV, \ \bth' \gets \bth$
$A_k$ in the operator pool
update $\bM''$ and $\bV''$ *[r]Add $A_k$ and $\theta_k\equiv 0$ to the ansätze
$\dot\bth'' \gets \bM'' \backslash \bV''$
$\mathcal{D}''\gets\dot{\bth''}^T \bM'' \dot{\bth''}-2 \bV'^T \dot{\bth''}$
$\mathcal{D}' \gets \mathcal{D}''$
$\bM'\gets\bM'',\ \bV'\gets\bV'', \ \bth'\gets \bth'', \ \dot{\bth}'\gets \dot{\bth}''$

$\bth\gets \bth',\ \mathbf{A}\gets \mathbf{A}',\ \dot{\bth}\gets \dot{\bth}'$

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 $\bth$ and $\mathbf{A}$ at each time step and keep track of the norm $e^{-\Gamma}$. The final vectorized state is $e^{-\Gamma/2}\ket{\phi\pqty{t_f}}$. However, to measure any observable $O$ in the unvectorized setting, we will need to evaluate

\begin{equation}
\expval{O} = e^{-\frac{\Gamma}{2}}\norm{O}_1\braket{O^\dagger}{\phi(t_f)}
\end{equation}

where $\ket{O^\dagger}$ is the normalized and vectorized operator $O^\dagger$ and $\norm{\cdot}_1$ denotes the trace norm.  This quantity can be measured either directly or indirectly with the help of a synthesized unitary $V$ which prepares $\ket{O^\dagger}$, i.e. $\ket{O^\dagger}=V\ket{0}^{\otimes 2N}$ <cit.> (See Appendix <ref> 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 <ref> 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 algorithm. It can be either a additive threshold or a multiplicative threshold. In this paper, we will use the additive threshold as described in Algorithm <ref>, and refer it as the adaptive threshold. In addition, the algorithms (restricted or unrestricted ) require us to solve the linear equations given by Eq. (<ref>), which can become ill-conditioned as the size of ansätze increases. To address this challenge, Tikhonov ($L_2$) regularization

\begin{equation}
\label{eq:tik_reg}
\dot{\bth} = \pqty{\bM^\dagger \bM + \lambda \mathbf{I}}^{-1}\bM^\dagger \bV \ ,
\end{equation}

is applied when inverting $\bM$. Here $\lambda$ is a small parameter to shift the diagonals of the $\bM^\dagger \bM$ matrix.

§ RESULT

§.§ 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
Ht = AtH_D + Bt H_P  ,

where $A\pqty{t}$ and $B\pqty{t}$ (known as the annealing schedule) are scalar functions which satisfy $A\pqty{0} \gg B\pqty{0} \approx 0$ and $B\pqty{t_f} \gg A\pqty{t_f} \approx 0$. Here, we will only consider the linear schedule
At = 1-t/t_f  ,   Bt = t / t_f  .

$H_\mathrm{D}$ and $H_\mathrm{P}$ are constant terms given by

\begin{align}
H_\mathrm{D} = -\sum_{i=1}^N X_i \ , \quad H_\mathrm{P} = -\sum_{i=1}^{N-1}j_iZ_iZ_{i+1} \ ,
\end{align}

where $N$ is the total number of qubits and $Z_i$ ($X_i$) is the Pauli $Z$ ($X$) matrix on the $i$th qubit. The couplings strength $j_i$s are given by

\begin{equation}
j_{i}= \begin{cases}w_{1} & \text { if }\lceil i / n\rceil \text { is odd } \\ w_{2} & \text { otherwise }\end{cases} \ ,
\end{equation}

where $n$ is the sector size. In this study, we fix the model parameters as $\omega_1 = 1$, $\omega_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:
L_i = Z_i  ,   i ∈1, N
with rate $\gamma_i$. It describes a continuous dephasing channel on each qubit. It is worth noting that the conditions outlined in Eqs. (<ref>) 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 $\gamma_i = 0.01$ in our simulations.

The second model consists of two Lindblad operators on each qubit
L_i^+ = ( X_i + i Y_i ) /2  ,  L_i^- = ( X_i - i Y_i ) /2
with different rate $\gamma_i^+$ and $\gamma_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 $\gamma_i^+ = 0.04$ and $\gamma_i^- = 0.004$ henceforth. It is worth noting that we choose the $\gamma$ 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 mid-circuit measurement or block encoding (unitary dilation) <cit.> based methods can be adopted to implement those non-unitary jumps (See Appendix <ref> for details).

We conducted numerical experiments to compare different operator pools.
The choice of our operator pool is a combination of the Hamiltonian pool <cit.> and the qubit-adapt pool proposed by Tang et. al. <cit.>. While there can be multiple choices of operator pools, variational ansätze generated with qubit-adapt pools are much shallower <cit.>. Starting from a pool consisting of all single qubit Pauli operators
$\mathcal{P}_{s} =\{X_i\}_{i=1}^{N}\cup\{Y_i\}_{i=1}^{N} \cup\{Z_i\}_{i=1}^{N}$, we constructed three distinct operator pools by adding three different types of two-qubit Pauli operators to $\mathcal{P}_{s}$, defined as follows:

\begin{align}
\mathcal{P}_{1} &= \mathcal{P}_{s}\cup \{Z_{i}Z_{i+1}\}_{i=1}^{N-1}\\
\mathcal{P}_{2} &= \mathcal{P}_{s}\cup \{P_{i}P_{i+1}\}_{i=1}^{N-1} \\
\mathcal{P}_{3} &= \mathcal{P}_{s}\cup \{P_{i}P_{j}\}_{i,j=1}^{N} \ ,
\end{align}

where $P_{i} \in \{ X_i, Y_i, Z_i\}_{i=1}^{N}$.  We use the individual terms from Eqs. (<ref>) as the operators in our pool.  $\mathcal{P}_{1}$ draws its inspiration from the Hamiltonian pool (the terms in the Hamiltonian appear in the pool only).
We found the $\mathcal{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 $\mathcal{P}_{2}$ generates the desired results for the trajectory method, it fails in the vectorization method case. Only $\mathcal{P}_{3}$ proves effective in producing the desired results for both methods. We summarize the pool dependence of our method in Table <ref>.

MethodPool     $\mathcal{P}_1$     $\mathcal{P}_2$     $\mathcal{P}_3$

Trajectory     $\times$     Yes     Yes

Vectorization     $\times$     $\times$     Yes

Dependence of our method on operator pools.

Noiseless simulation using the trajectory method and UAVQDS for  $\mathbf{N=8}$ ($\mathbf{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\sigma$ ($\sigma$ 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 $\mathcal{P}_2$.

Noiseless simulation using the vectorization method and UAVQDS for  $\mathbf{N=4}$ ($\mathbf{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. The adaptive threshold was set to $r=10^{-6}$ and the time step size was set to $dt=0.01$. We used the operator pool $\mathcal{P}_3$.

In Figs. <ref> and <ref>, 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 and eigenstate populations obtained from UAVQDS with the exact solution obtained from the Hamiltonian open quantum system toolkit (HOQST) package <cit.>. The eigenstate populations are chosen because they are the quantities of interest in the original QA setting, and also demonstrate that the UAVQDS correctly captures the dynamics of the high energy levels. While it is unknown whether the eigenstate populations can be efficiently measured on a real quantum computer, they can be trivially 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 $\lambda$ in Eq. (<ref>) to a fixed value of $10^{-8}$.

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, we utilized $r = 10^{-4}$ to generate the results for the $N=8$ problem (with $8$ physical qubits), while we had to decrease $r = 10^{-6}$ for the vectorization $N=4$ problem
(also with $8$ physical qubits) to ensure enough operators were included in the simulation (See Appendix <ref> 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. <ref>. To quantify the error, we define an infidelity measure $D$ as the trace-norm distance between the solution of UAVQDS and the exact solution at the end of the evolution, given by

\begin{equation}
D = \Tr \bqty{\sqrt{\pqty{\rho(t_f) - \rho_t}^\dagger\pqty{\rho(t_f)-\rho_t}}} \ ,
\end{equation}

where $\rho(t_f)$ and $\rho_t$ represent the solutions obtained from UAVQDS and the exact solver, respectively. As shown in Fig. <ref>, 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. <ref>. 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.

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 ($\mathcal{P}_2$). 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 ($\mathcal{P}_2$ and $\mathcal{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. <ref>.

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\theta P_l}$, where $P_l$ is a Pauli word of length $l$.  Thus,  the total number of CNOT gates needed is given by  $\sum_{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 Fig. <ref> and <ref> for both the trajectory and vectorization methods. 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 $a N^{b}$ or an exponential of the form $\alpha e^{\beta N}$ using least squares algorithm, and the results are summarized in Table. <ref>. 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 Fig. <ref> and <ref>.
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 <cit.>. Confirming this observation for larger system in case of open quantum systems will be a topic for future studies.

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. 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. <ref>.

DataFitting     $b$     $R^2$ ($\mathrm{poly}$)     $R^2$ ($\mathrm{exp}$)

Trajectory (median)     $2.39$     $0.986$     $0.961$

Trajectory (mean)     $1.98$     $0.980$     $0.938$

Trajectory (max)     $1.72$     $0.972$     $0.920$

Vectorization ($r=10^{-5}$)     $1.13$     $0.821$     $0.795$

Vectorization ($r=10^{-6}$)     $2.47$     $0.945$     $0.943$

A summary of fitting results. The ansatz counts shown in Fig. <ref> are fitted to either a polynomial of the form $a N^{b}$ or an exponential of the form $\alpha e^{\beta 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.

§.§ 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. (<ref>) on a real quantum computer. We run our experiments on IBM's 27-qubit processor  with the Falcon architecture. To ensure the algorithms run smoothly, it is necessary to carefully compile the circuit based on the specific quantum device being 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 <ref>). 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 <ref>). To synthesize the unitary required for measurements in the vectorization method, we use the Berkeley Quantum Synthesis Toolkit (BQSKit) <cit.>. It is worth noting that, at the time of this paper, the function to synthesize an arbitrary state preparation unitary was only available in the 1.1-dev branch of BQSKit. 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 <cit.> and the Qiskit documentation <cit.>. We also incorporate dynamical decoupling using Qiskit's standard tools, implementing periodic gate sequences to suppress cross-talk and system-environment coupling <cit.>. Additionally, we apply an empirical error mitigation scheme based on the classical resolution enhancement (RE) technique, which has shown effectiveness in a different context <cit.>, 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-the-box solutions and are resource-efficient, making them suitable for real-world applications.

Energy expectation value obtained by  executing the parameterized circuits on . 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 ($\mathrm{M3}$), dynamical decoupling (DD) and resolution enhancement (RE). The error bars in panel (a) represent two times the standard deviation from two 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.

We show the experimental results on the hardware in Fig. <ref>. The figure shows $\expval{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 $\expval{H}$ with the Hadamard test.

From Fig. <ref> (a) and (b)
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 $2$ 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. <ref>) 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 systems. Our algorithm pushes the boundary of NISQ algorithms for open quantum system simulation, providing a promising avenue for future research in this field.

The authors would like to thank Aaron Szasz for insightful discussions on BQSKit. This work was supported by the U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers, Quantum Systems Accelerator (HC, WAdJ) and by DOE Office of Advanced Scientific Computing Research (ASCR) through the ARQC program (NG).

This research used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725 and resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory.

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.

§ DERIVATION OF THE VECTORIZED LINDBLAD EQUATION

In this section we show how to derive the vectorized Lindblad equation from

\begin{equation}
\label{app:lind_eq}
\frac{d}{dt}\rho\pqty{t} = -i\comm{H\pqty{t}}{\rho\pqty{t}}+\Lag\bqty{\rho\pqty{t}} \ ,
\end{equation}

where the dissipative term $\Lag\bqty{\rho\pqty{t}}$ is

\begin{equation}
\Lag\bqty{\rho\pqty{t}} = \sum_{k=1}^{K}  \pqty{L_k \rho\pqty{t} L^\dagger_k -\frac{1}{2}\acomm{L^\dagger_k L_k}{\rho\pqty{t}}} \ .
\end{equation}

The strategy is to apply the identity $\operatorname{vec}(A B C)=(C^{\mathrm{T}} \otimes A)\operatorname{vec}(B)$ to each term on the right-hand side of Eq. (<ref>). The Hamiltonian part $\comm{H\pqty{t}}{\rho\pqty{t}}$ can be written in terms of

\begin{equation}
\label{app:Hvec}
\operatorname{vec} \Bqty{  \rho\pqty{t} H\pqty{t}} = \pqty{H^\mathrm{T}\pqty{t}\otimes I}\operatorname{vec}\pqty{\rho\pqty{t}}, \quad
\operatorname{vec} \Bqty{ H\pqty{t} \rho\pqty{t}} = \pqty{I\otimes H\pqty{t}}\operatorname{vec}\pqty{\rho\pqty{t}} \ .
\end{equation}

The Lindblad part can be written in terms of

\begin{equation}
\label{app:Lvec}
\operatorname{vec} \pqty{L_k \rho L^\dagger_k} = \pqty{L_k^*\otimes L_k }\operatorname{vec}\pqty{\rho} ,\quad \operatorname{vec} \pqty{L^\dagger_k L_k \rho } = \pqty{I \otimes L^\dagger_k L_k }\operatorname{vec}\pqty{\rho} ,\quad \pqty{ \rho L^\dagger_k L_k} =  \pqty{L^\mathrm{T}_k L^*_k \otimes I} \operatorname{vec}\pqty{\rho} \ .
\end{equation}

Summing up every term in Eqs.(<ref>) and (<ref>), the effective Hamiltonian can be obtained as

\begin{align}
H_\mathrm{eff} = I\otimes H - H^T\otimes I
+i\sum_k \bqty{L_k^*\otimes L_k-\frac{1}{2}\pqty{I\otimes L_k^{\dagger} L_k + L_k^{T} L^*_k\otimes I}} \ .
\end{align}

§ DERIVATION OF THE EVOLUTION EQUATION

In this appendix, we derive the evolution equation of the variational parameters for the effective Schrödinger equation

|*⟩ψ̇ = -i H_eff|ψ⟩  ,
where the dot symbol denotes the time derivative. Before proceeding, it is worth mentioning that the effective Hamiltonian in Eq. (<ref>) is not necessarily Hermitian. The ansatz state is parameterized by $\bth\pqty{t}\equiv\bqty{\theta_1\pqty{t}, \theta_2\pqty{t}, \cdots \theta_k\pqty{t}}^T$

\begin{equation}
\label{app_eq:ansatz}
\ket{\phi\pqty{\bth\pqty{t}}} = \prod_{\mu=1}^{k} e^{-i\theta_\mu\pqty{t} A_\mu} \ket{\psi_\mathrm{R}} \ ,
\end{equation}

where $A_\mu$ and $\ket{\psi_\mathrm{R}}$ are the ansatz operators and the reference state, respectively.
It is more convenient to derive the evolution equation of $\bth\pqty{t}$ using the density matrix representation to avoid the problem of a time dependent global phase <cit.>. Let $W\pqty{\bth}\equiv\dyad{\phi\pqty{\bth}}{\phi\pqty{\bth}}$, the equation of motion is given by

\begin{align}
\dot{W} &= \dyad*{\dot{\phi}}{\phi}+\dyad*{\phi}{\dot{\phi}} \label{app_eq:Wdot}\\
&= -i\pqty{\Heff W - W \Heff^\dagger} \label{app_eq:Wcomm}\ ,
\end{align}

where we make use of Eq. (<ref>)
in going from Eq. (<ref>) to Eq. (<ref>). Applying the McLachlan’s variational principle to Eq. (<ref>), we have

δdWt/dt-Wt^2 =0  ,
where $\Lag\bqty{W}\equiv  -i\pqty{\Heff W - W \Heff^\dagger}$.
Let $\mathcal{D}\equiv\norm{\frac{dW\pqty{\bth\pqty{t}}}{dt}-\Lag\bqty{W\pqty{\bth\pqty{t}}}}^{2}$, the McLachlan distance (which is actually $\sqrt{\mathcal{D}}$) can be simplified as

𝒟=Ẇ^†-^†WẆ-W = Ẇ^†Ẇ -2ReẆ^†W+^†WW  .

Based on Eq. (<ref>), $\dot{W}$ is Hermitian, i.e., $\dot{W}=\dot{W}^\dagger$. The first term $\Tr\bqty{\dot{W}^\dagger \dot{W}}$ on the RHS of Eq. (<ref>) can be simplified by the chain rule
Ẇ^†Ẇ = ∑_μ,νWθ_μWθ_ν θ̇_μθ̇_ν= ^T   ,
where $\bM$ is a matrix with elements

\begin{align}
\bM_{\mu \nu} &=\Tr\Bqty{\pdv{W}{\theta_\mu}\pdv{W}{\theta_\nu}} \\
&=\Tr\Bqty{\pqty{ \pdv{\ket{\phi}}{\theta_\mu}\bra{\phi} + \ket{\phi}\pdv{\bra{\phi}}{\theta_\mu} }\pqty{ \pdv{\ket{\phi}}{\theta_\nu}\bra{\phi} + \ket{\phi}\pdv{\bra{\phi}}{\theta_\nu} }}\\
&= 2\mathrm{Re}\pqty{\pdv{\bra{\phi}}{\theta_{\mu}}\pdv{\ket{\phi}}{\theta_{\nu}} +\bra{\phi}\pdv{\ket{\phi}}{\theta_{\mu}} \bra{\phi}\pdv{\ket{\phi}}{\theta_{\nu}}} \ .
\end{align}

Because $\theta_\mu$ are real parameters, $\bM$ is a real symmetric matrix.

Similarly, the second term on the RHS of Eq. (<ref>) can be simplified as

\begin{align}
2\mathrm{Re}\Bqty{\Tr\Bqty{\dot{W}^\dagger\Lag\bqty{W}}} &= 2\mathrm{Re}\Bqty{\Tr\Bqty{\pdv{W}{\theta_\mu} \Lag\bqty{W}}}\dot{\theta}_\mu \\
&=2 \bV^T \dot{\bth} \ ,
\end{align}

where $\bV$ is a vector with elements

\begin{align}
\bV_\mu &= \mathrm{Re}\Bqty{\Tr\Bqty{\pdv{W}{\theta_\mu}\Lag\bqty{W}}}\\
&=\mathrm{Im}\Bqty{\Tr\Bqty{ \pqty{ \pdv{\ket{\phi}}{\theta_\mu}\bra{\phi} + \ket{\phi}\pdv{\bra{\phi}}{\theta_\mu}} \pqty{\Heff \dyad{\phi}{\phi}-\dyad{\phi}{\phi} \Heff^\dagger}}} \\
&=2\mathrm{Im}\pqty{\expval{\Heff}\bra{\phi}\pdv{\ket{\phi}}{\theta_\mu}+\pdv{\bra{\phi}}{\theta_\mu}\Heff \ket{\phi}} \ ,
\end{align}

and $\expval{\Heff} = \mel{\phi}{\Heff}{\phi}$.

Before calculating the last term on the RHS of Eq. (<ref>), we split the effective Hamiltonian into a Hermitian and an anti-Hermitian parts

H_eff = H_eff+H^†_eff/2 -i iH_eff-H^†_eff/2 ≡H_e - i H_a  .
Then the last term becomes

\begin{align}
\Tr\Bqty{\Lag^\dagger\bqty{W}\Lag\bqty{W}} &= \Tr\Bqty{\pqty{W\Heff^\dagger-\Heff W}\pqty{\Heff W - W \Heff^\dagger}} \\
&=\Tr\bqty{W\Heff^\dagger\Heff W}-\Tr\bqty{W\Heff^\dagger W \Heff^\dagger} - \Tr\bqty{\Heff W \Heff W} + \Tr\bqty{\Heff W W \Heff^\dagger} \\
&=\expval{\pqty{\He + i \Ha}\pqty{\He -i\Ha}} - \expval{\He + i\Ha}^2 - \expval{\He - i\Ha}^2 + \expval{\pqty{\He+i\Ha}\pqty{\He-i\Ha}} \\
&=2\expval{\He^2} - 2\expval{\He}^2 + 2\expval{\Ha^2} + 2\expval{\Ha}^2 +2i\expval{\comm{\Ha}{\He}} \ .
\end{align}

The variational principle (Eq. (<ref>)) then yields
δ𝒟 = ^T + ^T - 2^T = 2 - 2=  0  .
Finally, the evolution equation of the variational parameters is given by a linear equation
=  .

§ LOWER BOUND OF THE MCLACHLAN'S DISTANCE

In this Appendix, we derive a lower bound of the McLachlan's distance (Eq. (<ref>)). First, we note from Eq. (<ref>) that the ansatz state is parameterized by a unitary

\begin{equation}
\ket{\phi\pqty{\bth\pqty{t}}} = U\pqty{\bth\pqty{t}} \ket{\psi_\mathrm{R}} \ .
\end{equation}

A reasonable assumption we make is that $d U\pqty{\bth\pqty{t}} / dt$ exists. With this assumption, we can define an effective Hamiltonian parameterized by the derivative of the ansatz parameter
H̃ = i U̇ U^-1  ,
where we omit the $t$-dependence of each quantity in the above equation. It is straightforward to check the following relation is satisfied

U̇ = -i H̃ U  .
Then the derivative of the ansatz state is subject to the following equation of motion
dW/dt = - iH̃W  ,
and the McLachlan's distance becomes

\begin{align}
\mathcal{D} &= \norm{\frac{dW\pqty{\bth\pqty{t}}}{dt}-\Lag\bqty{W\pqty{\bth\pqty{t}}}}^{2}    \\
&=\norm{\pqty{\tilde{H}-\Heff}W-W\pqty{\tilde{H}-\Heff^\dagger}}^2 \label{app_eq:H-Heff}\\
&=\norm{\pqty{\Delta H+i\Ha}W-W\pqty{\Delta H-i\Ha}}^2 \label{app_eq:DH-Ha}\\
&=2\expval{\Delta H^2} - 2\expval{\Delta H}^2 + 2\expval{\Ha^2} + 2\expval{\Ha}^2 -2i\expval{\comm{\Ha}{\Delta H}} \ ,
\end{align}

where we denote $\Delta H\equiv \tilde{H}-\He$ in going from Eq. (<ref>) to Eq. (<ref>).
When $\comm{\Ha}{\Delta H}=0$, a lower bound of $\mathcal{D}$ is given by
𝒟 ≥2^2 + 2^2  ,
where the equality is achieved when $\Delta H=0$, which means the optimal $\dot{\bth}$ makes $\tilde{H}\pqty{\dot{\bth}}$ in Eq. (<ref>) equal the Hermitian part of the effective Hamiltonian $\Heff$.

A sufficient condition for $\comm{\Ha}{\Delta H}=0$ is

\begin{align}
\comm{\He}{\Ha} &= 0 \\
\comm{A_\mu}{\Ha} & = 0 \ , \forall \mu \in \bqty{1, k} \label{app_eq:comm_A}\ ,
\end{align}

where Eq. (<ref>) guarantees $\comm{\tilde{H}\pqty{\dot{\bth}}}{\Ha}=0$ .

§ 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

\begin{equation}
\frac{d}{dt}\ket{\tilde{\psi}\pqty{t}} = -i H_\mathrm{eff}\ket{\tilde{\psi}\pqty{t}} \ ,
\end{equation}

with $\ket{\tilde{\psi}\pqty{t}}$ being the unnormalized state vector and $\ket{\psi\pqty{t}}$ being the normalized state vector. The effective Hamiltonian has both the Hermitian $H_\mathrm{e}$ and anti-Hermitian parts $-iH_\mathrm{a}$, i.e.,
H_\mathrm{eff} = H_\mathrm{e} - i H_\mathrm{a}
$. We also use the notation $\ket{\phi\pqty{t}}$ for the state generated by our variational circuit. The evolution of $\braket{\tilde{\psi}\pqty{t}}$ can be calculated as
\begin{subequations}
\begin{align}
\frac{d\braket{\tilde{\psi}\pqty{t}}}{dt} &= \frac{d\bra{\tilde\psi \pqty{t}}}{dt} \ket{\tilde\psi\pqty{t}} + \bra{\tilde\psi\pqty{t}}\frac{d\ket{\tilde\psi\pqty{t}}}{dt} \\
&= \mel{\tilde\psi\pqty{t}}{\pqty{iH^\dagger_\mathrm{eff} - iH_\mathrm{eff}}}{\tilde\psi\pqty{t}}\\
&= -2 \mel{\tilde\psi\pqty{t}}{H_\mathrm{a}}{\tilde\psi\pqty{t}}\\
&= -2 \mel{\psi\pqty{t}}{H_\mathrm{a}}{\psi\pqty{t}} \braket{\tilde\psi\pqty{t}} \ .
\end{align}
\end{subequations}
The solution to the above equation is
\begin{equation}
\braket{\tilde\psi\pqty{t}} = e^{-2 \int_0^{t} \mel{\psi\pqty{\tau}}{H_\mathrm{a}}{\psi\pqty{\tau}} d\tau} \braket{\tilde\psi\pqty{0}} \ .
\end{equation}
In the above equation, $e^-2 ∫_0^t ψτH_aψτ dτ$ is monotonically decreasing. To see this, we first examine the case of the unravelled Lindblad equation, where $H_a = 1/2∑_k L^†_k L_k$ is positive semi-definite since each term in the summation is quadratic. As a result, $∫_0^t ψτ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.,
\begin{equation}
\sqrt{\Tr\bqty{\rho^\dagger\pqty{t_1}\rho\pqty{t_1}}} \leq \sqrt{\Tr\bqty{\rho^\dagger\pqty{t_2}\rho\pqty{t_2}}} \ , \quad \mathrm{for} \ t_1 \geq t_2 \ ,
\end{equation}
and the linear isometry preserves the trace norm and $L_2$ norm, we have
\begin{equation}
\braket{\tilde\psi\pqty{t_1}} \leq \braket{\tilde\psi\pqty{t_2}} \ , \quad \mathrm{for} \ t_1 \geq t_2 \ .
\end{equation}

For a small time interval $dt$
\begin{equation}
\braket{\tilde\psi\pqty{t+dt}} \approx e^{-2 \mel{\psi\pqty{t}}{H_\mathrm{a}}{\psi\pqty{t}} dt} \braket{\tilde\psi\pqty{t}} \ ,
\end{equation}
thus we can keep track of the norm by measuring $H_a$ at each time step.
\section{Measurement circuit for $\bM$ and $\bV$}
\label{app:MV_measurement}
In this appendix, we will discuss the quantum circuit used to measure the parameter in the variational equation of motion $ = $, where
\be
\label{app:M}
\bM_{\mu \nu} =2\mathrm{Re}\pqty{\pdv{\bra{\phi}}{\theta_{\mu}}\pdv{\ket{\phi}}{\theta_{\nu}} +\bra{\phi}\pdv{\ket{\phi}}{\theta_{\mu}} \bra{\phi}\pdv{\ket{\phi}}{\theta_{\nu}}} \ ,
\ee
\be
\label{app:V}
\bV_\mu = 2\mathrm{Im}\pqty{\expval{\Heff}\bra{\phi}\pdv{\ket{\phi}}{\theta_\mu}+\pdv{\bra{\phi}}{\theta_\mu}\Heff \ket{\phi}} \ .
\ee
assuming that each ansatz operator $A_μ$ is a Pauli operator. We only present the resulting circuits here and encourage interested readers to refer to~\cite{Mitarai2019-pv} 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.~\ref{fig:measure_pauli}, where we use notation $Uθ⃗$ to denote the the unitary generated by the variational circuit and $HS^†$ to denote the optional Hadamard (Hadamard-phase) gate to rotate the basis of $P_i$ to $z$ direction.
\begin{figure}[htb]
\centering
\subfloat[]{
\begin{quantikz}
\lstick{$\ket{0}$} & \gate{H} & \ctrl{1} & \gate{H} & \meter{} \\
\lstick{$\ket{\psi_\mathrm{R}}$} &  \gate{U(\vec\theta)}\qwbundle[
alternate]{} & \gate{P_i} \qwbundle[
alternate]{} & \qwbundle[
alternate]{} & \qwbundle[
\end{quantikz}
\qquad
\subfloat[]{%
\begin{quantikz}
\lstick{$\ket{\psi_\mathrm{R}}$} &  \gate{U(\vec\theta)}\qwbundle[
alternate]{} & \gate{H(S^\dagger)} \qwbundle[
\gategroup[1,steps=2, style={dashed,
rounded corners, fill=blue!20, inner xsep=2pt}, background]{{$\mathcal{M}_{P_i}$}}& \meter{}\qwbundle[
\end{quantikz}
\caption{Quantum circuit for evaluating $\ev{P_i}{\phi}$ using (a) Hadamard test; (b) direct measurements. $U\pqty{\vec\theta}$ denotes the unitary generated by the variational circuit $\prod_{\mu=1}^{k} e^{-i\theta_\mu\pqty{t} A_\mu}$. $H\pqty{S^\dagger}$ 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$.}
\label{fig:measure_pauli}
\end{figure}

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
\begin{equation}
\pdv{\ket{\phi}}{\theta_\mu} = -i \prod_{l=\mu+1}^{k} e^{-i\theta_l A_l}  A_\mu \prod_{j=1}^{\mu} e^{-i\theta_j A_j} \ket{\psi_\mathrm{R}} \ ,
\end{equation}
the term $⟨ϕ|P_i|ϕ⟩θ_μ$ can be simplified as
\begin{equation}
\label{app:PA_measure}
\bra{\phi}P_i\pdv{\ket{\phi}}{\theta_\mu} \propto \ev{U^\dagger_{\mu:1}U^\dagger_{k:\mu+1}P_i U_{k:\mu+1} A_\mu U_{\mu:1}}{\psi_\mathrm{R}} \ ,
\end{equation}
where $U_μ:1 = ∏_j=1^μ e^-iθ_j A_j$ and $U_k:μ+1 = ∏_l=μ+1^k e^-iθ_l A_l$. Eq.~\eqref{app:PA_measure} can be evaluated using the generalized Hadamard test (Fig.~\ref{fig:g_hadamard}) or the corresponding direct measure circuit (Fig.~\ref{fig:g_dm1} and~\ref{fig:g_dm2}).

\begin{figure}[htb]
\centering
\subfloat[]{
\begin{quantikz}
\lstick{$\ket{0}$} & \gate{H} & \ctrl{1} & \qw & \ctrl{1} & \gate{S^b} & \gate{H} & \meter{} \\
\lstick{$\ket{\psi_\mathrm{R}}$} &  \gate{U_{\mu:1}}\qwbundle[
alternate]{} & \gate{A_\mu} \qwbundle[
alternate]{} & \gate{U_{k:\mu+1}}\qwbundle[
alternate]{} &\gate{P_i}\qwbundle[
alternate]{} &\qwbundle[
alternate]{} &\qwbundle[
alternate]{} &\qwbundle[
\end{quantikz}
\label{fig:g_hadamard}
\\
\subfloat[]{%
\begin{quantikz}
\lstick{$\ket{\psi_\mathrm{R}}$} &  \gate{U_{\mu:1}}\qwbundle[
alternate]{} & \gate{H(S^\dagger)} \qwbundle[
\gategroup[1,steps=2, style={dashed,
rounded corners, fill=blue!20, inner xsep=2pt}, background]{{$\mathcal{M}_{A_\mu}$}}& \meter{}\qwbundle[
alternate]{}& \gate{U_{k:\mu+1}}\qwbundle[
alternate]{} & \gate{H(S^\dagger)} \qwbundle[
\gategroup[1,steps=2, style={dashed,
rounded corners, fill=blue!20, inner xsep=2pt}, background]{{$\mathcal{M}_{P_i}$}}& \meter{}\qwbundle[
\end{quantikz}
\label{fig:g_dm1}
\\
\subfloat[]{
\begin{quantikz}
\lstick{$\ket{\psi_\mathrm{R}}$} &  \gate{U_{\mu:1}}\qwbundle[
alternate]{} & \gate{e^{\pm i \pi A_\mu / 4}} \qwbundle[
& \gate{U_{k:\mu+1}}\qwbundle[
alternate]{} & \gate{H(S^\dagger)} \qwbundle[
\gategroup[1,steps=2, style={dashed,
rounded corners, fill=blue!20, inner xsep=2pt}, background]{{$\mathcal{M}_{P_i}$}}& \meter{}\qwbundle[
\end{quantikz}
\label{fig:g_dm2}
\caption{Quantum circuits for evaluating $\bra{\phi}P_i\pdv{\ket{\phi}}{\theta_\mu}$ using (a) generalized Hadamard test; (b) and (c) direct measurements. The generalized Hadamard test circuit uses a binary integer $b \in \{0,1\}$ as an input. If $b=0$, the circuit measures the real part of $\ev{U^\dagger_{\mu:1}U^\dagger_{k:\mu+1}P_i U_{k:\mu+1} A_\mu U_{\mu:1}}{\psi_\mathrm{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.}
\label{fig:measure_gpauli}
\end{figure}

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
\begin{equation}
\pdv{\bra{\phi}}{\theta_{\mu}}\pdv{\ket{\phi}}{\theta_{\nu}} \propto \ev{U^\dagger_{\mu:1}A_\mu U_{\mu+1:\nu+1} A_\nu U_{\nu:1}}{\psi_\mathrm{R}} \ ,
\end{equation}
assuming $μ≥ν+1$. Quantum circuits for evaluating this quantity are shown in Fig.~\ref{fig:measure_metric_tensor}.
\begin{figure}[htb]
\centering
\subfloat[]{
\begin{quantikz}
\lstick{$\ket{0}$} & \gate{H} & \ctrl{1} & \gate{X}\qw & \ctrl{1} & \gate{X}\qw & \gate{S^b} & \gate{H} & \meter{} \\
\lstick{$\ket{\psi_\mathrm{R}}$} &  \gate{U_{\nu:1}}\qwbundle[
alternate]{} & \gate{A_\nu} \qwbundle[
alternate]{} & \gate{U_{\mu:\nu+1}}\qwbundle[
alternate]{} &\gate{A_\mu}\qwbundle[
alternate]{} &\qwbundle[
alternate]{} &\qwbundle[
alternate]{} &\qwbundle[
alternate]{} &\qwbundle[
\end{quantikz}
\\
\subfloat[]{%
\begin{quantikz}
\lstick{$\ket{\psi_\mathrm{R}}$} &  \gate{U_{\nu:1}}\qwbundle[
alternate]{} & \gate{H(S^\dagger)} \qwbundle[
\gategroup[1,steps=2, style={dashed,
rounded corners, fill=blue!20, inner xsep=2pt}, background]{{$\mathcal{M}_{A_\nu}$}}& \meter{}\qwbundle[
alternate]{}& \gate{U_{\mu:\nu+1}}\qwbundle[
alternate]{} & \gate{H(S^\dagger)} \qwbundle[
\gategroup[1,steps=2, style={dashed,
rounded corners, fill=blue!20, inner xsep=2pt}, background]{{$\mathcal{M}_{A_\mu}$}}& \meter{}\qwbundle[
\end{quantikz}
\\
\subfloat[]{
\begin{quantikz}
\lstick{$\ket{\psi_\mathrm{R}}$} &  \gate{U_{\nu:1}}\qwbundle[
alternate]{} & \gate{e^{\pm i \pi A_\nu / 4}} \qwbundle[
& \gate{U_{\mu:\nu+1}}\qwbundle[
alternate]{} & \gate{H(S^\dagger)} \qwbundle[
\gategroup[1,steps=2, style={dashed,
rounded corners, fill=blue!20, inner xsep=2pt}, background]{{$\mathcal{M}_{A_\mu}$}}& \meter{}\qwbundle[
\end{quantikz}
\caption{Quantum circuits for evaluating $\pdv{\bra{\phi}}{\theta_{\mu}}\pdv{\ket{\phi}}{\theta_{\nu}}$ using (a) generalized Hadamard test; (b) and (c) direct measurements. The generalized Hadamard test circuit uses a binary integer $b \in \{0,1\}$ as an input. If $b=0$, the circuit measures the real part of $\ev{U^\dagger_{\mu:1}A_\mu U_{\mu+1:\nu+1} A_\nu U_{\nu:1}}{\psi_\mathrm{R}}$ (for $\mu \ge \nu$), 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.}
\label{fig:measure_metric_tensor}
\end{figure}
\section{Measure observable $O$ in the vectorization method}
\label{app:mO}
In this section we describe the circuit used to measure the observable $O$ for the vectorization method based on Ref.~\cite{Kamakari2022-yr}. 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
\begin{equation}
\ket{\phi\pqty{t}} = \prod_{\mu=1}^{k} e^{-i\theta_\mu\pqty{t} A_\mu} \ket{\psi_\mathrm{R}} = U\pqty{\bth} \ket{\psi_\mathrm{R}} =U\pqty{\bth} U_\mathrm{R}\ket{0} \ ,
\end{equation}
where $U_R$ prepares the reference state.
If we can construct another unitary that prepares $|O^†⟩$, i.e. $|O^†⟩ = V |0⟩$, then $⟨O^†|ϕ(t)⟩=0V^†UU_R0$ 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
\begin{equation}
\label{app:vec_O}
\ket{O^\dagger}= \sum_{x_1 y_1}\frac{O^*_{x_1 y_1}}{\sqrt{\Tr\pqty{O^\dagger O}}}\ket{x_1y_1} \ ,
\end{equation}
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~\cite{Sun2023-mc,szasz2023numerical}. Because $O$ is usually a local observable, Eq.~\eqref{app:vec_O} can be further simplified to avoid any exponential cost~\cite{Kamakari2022-yr}.  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)~\cite{Younis2021-qs,szasz2023numerical} to synthesize the unitary $V$ and used the direct measurement protocol shown in Fig.~\ref{fig:vectorized_measure} 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)⟩ = V^†Uθ⃗ H^⊗2N0$. 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.
\begin{figure}[htb]
\begin{quantikz}
\lstick{$\ket{0}$} &  \gate{H^{\otimes 2N}}\qwbundle[
alternate]{} & \gate{U(\vec{\theta})} \qwbundle[
& \gate{V^\dagger}\qwbundle[
alternate]{} & \meter{}\qwbundle[
\end{quantikz}
\caption{Direct measurement scheme to evaluate $\braket{O^\dagger}{\phi(t)}$. The final result can be obtained by taking the square root of the all-$0$ string frequency.}
\label{fig:vectorized_measure}
\end{figure}
\section{Implementing $S_-$ and $S_+$ jump operators}
\label{app:jumps}
In this appendix, we describe two methods for applying the $S_+$ and $S_-$ jump operators on a quantum computer. Without loss of generality, we present only the derivation for $S_+≡(X+iY)/2$, and the procedure is easily extended to $S_-≡(X-iY)/2$. The key observation is that the $S_+$ jump only occurs when its probability $p_+∝S_-S_+ϕ$ is greater than zero, meaning $|ϕ⟩≠|0⟩$. Under this condition, the action of the jump operator $S_+|ϕ⟩/√(S_-S_+ϕ)$ is equivalent to a quantum eraser channel $ℰρ=0$. There are two ways of implementing such a channel on a quantum computer: mid-circuit measurement or block encoding (unitary dilation)~\cite{Low2017-ye,camps2023explicit,Del_Re2020-ul}. 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 Fig.~\ref{fig:q_eraser}.
\begin{figure}[htb]
\centering
\subfloat[]{
\begin{quantikz}
\lstick{$\ket{\phi}$} & \meter{$I/X$} \arrow[r] & \rstick{$\ket{0}$}
\end{quantikz}
\qquad
\subfloat[]{%
\begin{quantikz}
\lstick{$\ket{0}$}& \targ{} & \meter{} \arrow[r]& \rstick{$\ket{0}$} \\
\lstick{$\ket{\phi}$}& \ctrl{-1}  &\qw &\qw
\end{quantikz}
\caption{Quantum circuit for implementing an eraser channel $\mathcal{E}\pqty{\rho}=\dyad{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 $\ket{0}$ in (b) represents the post-selection. We only keep the runs where $0$ is returned by the measurement.}
\label{fig:q_eraser}
\end{figure}
\section{Growth of ansatz size}
\label{app:growth_ansatz}
In this appendix, we demonstrate how the sizes of ans\"atze grow as we adaptively add operators to them. In Fig.\ref{fig:traj_N_G}, 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\"atze after they are reset by the jumps. Consequently, the accumulated ansatz size becomes larger when more jumps occur. In Fig.\ref{fig:vec_N_G}, we showcase examples using the vectorization method with different error thresholds. We observe that larger ansatz sizes when lower threshold values are used.
\begin{figure}[htb]
\centering
\subfloat[]{
\label{fig:traj_N_G}
\includegraphics[width=0.47\linewidth]{figure/ASC_8_qubit_lind_Z_rcut_1.00e-04_unrestricted_parameters_num.pdf}
\hfill
\subfloat[]{%
\label{fig:vec_N_G}
\includegraphics[width=0.47\linewidth]{figure/ASC_4_qubit_lind_Z_vectorized_unrestricted_parameters_num.pdf}
\caption{\textbf{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\"atze 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.}
\label{fig:n4_param_grow}
\end{figure}

\section{The vectorization method with different operator pools and error threshold}
\label{app:vec_figs}
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 visual representations of how the error behaves with respect to various operator pools and error thresholds.
\begin{figure}[h]
\centering
\subfloat[]{
\includegraphics[width=0.45\linewidth]{figure/ASC_4_qubit_lind_Z_vectorized_unrestricted_err_H.pdf}
\hfill
\subfloat[]{%
\includegraphics[width=0.45\linewidth]{figure/ASC_4_qubit_lind_PM_vectorized_unrestricted_all_err_H.pdf}
\\
\subfloat[]{
\includegraphics[width=0.45\linewidth]{figure/ASC_4_qubit_lind_Z_vectorized_unrestricted_neighbor_err_H.pdf}
\hfill
\subfloat[]{
\includegraphics[width=0.45\linewidth]{figure/ASC_4_qubit_lind_PM_vectorized_unrestricted_err_H.pdf}
\caption{\textbf{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 ($\mathcal{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 ($\mathcal{P}_2$ in the main text) for the dephasing and amplitude damping models, respectively.}
\label{fig:n4_err_demo}
\end{figure}

\end{document}