[1]Institute of Physics, Ecole Polytechnique Fédérale de Lausanne (EPFL),
CH-1015 Lausanne, Switzerland [2]Center for Quantum Science and Engineering,
Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland

# Low-Rank Variational Quantum Algorithm for the Dynamics of Open Quantum
Systems

Sara Santos<EMAIL_ADDRESS>Xinyu Song<EMAIL_ADDRESS>Vincenzo Savona<EMAIL_ADDRESS>
###### Abstract

The simulation of many-body open quantum systems is key to solving numerous
outstanding problems in physics, chemistry, material science, and in the
development of quantum technologies. Near-term quantum computers may bring
considerable advantage for the efficient simulation of their static and
dynamical properties, thanks to hybrid quantum-classical variational
algorithms to approximate the dynamics of the density matrix describing the
quantum state in terms of an ensemble average. Here, a variational quantum
algorithm is developed to simulate the real-time evolution of the density
matrix governed by the Lindblad master equation, under the assumption that the
quantum state has a bounded entropy along the dynamics, entailing a low-rank
representation of its density matrix. The algorithm encodes each pure state of
the statistical mixture as a parametrized quantum circuit, and the associated
probabilities as additional variational parameters stored classically, thereby
requiring a significantly lower number of qubits than algorithms where the
full density matrix is encoded in the quantum memory. Two variational Ansätze
are proposed, and their effectiveness is assessed in the simulation of the
dynamics of a 2D dissipative transverse field Ising model. The results
underscore the algorithm’s efficiency in simulating the dynamics of open
quantum systems in the low-rank regime with limited quantum resources on a
near-term quantum device.

##  1 Introduction

Quantum systems are invariably influenced by their environment [1], an
interaction that often has detrimental effects such as decoherence, thereby
limiting the efficacy of quantum technological platforms. The ability to
efficiently simulate open quantum systems is thus of paramount importance,
providing essential insights into their fundamental properties and aiding in
the development of quantum technologies.

The statistical properties of an open quantum system are fully described by
the density matrix, which gives direct access to statistical ensemble averages
of physical quantities. The time evolution of the density matrix is governed
by a master equation which, within the general assumption of a Markovian
environment, takes the celebrated Lindblad form [1]. Alternatively, open
quantum systems admit a stochastic description in terms of quantum
trajectories, which simulate their dynamics by unraveling the evolution of the
system’s density matrix into an ensemble of stochastically evolving pure
states [2, 3, 4, 5]. A variety of numerical methods have been developed for
the simulation of open quantum system dynamics, either through direct
integration of the master equation or by averaging over simulated quantum
trajectories [6]. One of the main challenges in the simulation of open quantum
systems is the increased computational complexity due to the mixed nature of
the quantum state. If the unitary dynamics of an $n$-qubit system is described
by the time-evolution of its $N=2^{n}$ amplitudes, a mixed state of the same
system requires expressing $N^{2}$ elements of the density matrix. Simulating
one quantum trajectory has a computational complexity comparable to pure
quantum states, but averaging over many trajectories is required to estimate
physical quantities with sufficient accuracy. Finally, it has been shown [7]
that the entanglement entropy of the density matrix grows in time
significantly more slowly than that of individual quantum trajectories,
thereby making the approximate numerical description more efficient.

Quantum computers excel at representing quantum states efficiently, offering
the potential to simulate the real-time evolution of quantum systems, as
initially proposed by Richard Feynman [8]. Seth Lloyd has shown that the time
evolution operator $\hat{U}(t)=e^{-i\hat{H}t}$ can be implemented on quantum
computers efficiently when the Hamiltonian $\hat{H}$ of the system only
carries local interaction [9]. Quantum algorithms based on the Suzuki–Trotter
decomposition of the operator $\hat{U}(t)$ have become an active research
field [10, 11]. However, near-term, a.k.a NISQ devices [12], with limited
coherence time and operation accuracy, put constraints on the application of
Suzuki-Trotter time-evolution schemes, since the required quantum circuit
depth scales proportionally to the time length. To address the limitations of
NISQ devices, variational quantum algorithms (VQAs), have emerged [13, 14,
15]. VQAs efficiently prepare quantum states as low-depth parameterized
quantum circuits and evolve/optimize the parameters to capture the non-
equilibrium dynamics of a given physical system.

Several quantum algorithms for simulating the dynamics of open quantum systems
have been proposed [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29,
30, 31, 32, 33, 34]. These are mostly based on an appropriate dilation of the
system’s Hilbert space to encode environment degrees of freedom, and on mid-
circuit measurements and/or state reinitializations to model dissipative
processes, or rely on the Suzuki-Trotter expansion of the non-unitary dynamics
[29, 30]. Recently, variational quantum algorithms have been introduced to
simulate the Lindblad master equation directly, relying on a full encoding of
the density matrix as a parametrized circuit on a $2n$-qubit quantum register
[27, 31, 33, 34]. Similar parameterized Ansätze for learning mixed and thermal
states [35, 36, 37, 38] have also been studied.

Under several common conditions, open quantum systems are characterized by low
entropy, corresponding to quantum states with a low degree of mixedness.
Typical cases are the transient dynamics starting from an initial pure state,
weak system-bath coupling, as well as many quantum gate protocols on quantum
information platforms such as bosonic quantum codes, among others. A low-
entropy quantum state admits a low-rank approximation, by retaining only the
dominant density matrix eigenvalues. The low-rank approximation is therefore a
valuable opportunity to significantly reduce the computational complexity of
the simulation. Several low-rank approaches have been recently developed for
the simulation on classical computers [39, 36, 40], including a time-dependent
variational algorithm with adaptive dynamical rank [41].

Here we develop a variational quantum algorithm that simulates the dynamics of
open quantum systems with a truncated rank parameterized density matrix. The
mixed state is represented as a statistical mixture of pure states, and each
pure state is encoded as a parametrized quantum circuit, while the
probabilities, also entering as variational parameters, are stored
classically, taking advantage of the low-rank assumption of
$O(\mathrm{poly(n)})$ nonzero probabilities. The main advantage of this
representation is that the required quantum resources are the same as for a
pure quantum state, e.g. $n$ qubits for a system of $n$ interacting spins, as
opposed to the $2n$ qubits that would be required to encode the most general
density matrix. We introduce two variational Ansätze, characterized by
different trade-offs between expressive power and required computational
resources.

We demonstrate the algorithm on a 2D dissipative transverse field Ising model
(TFIM) and discuss the practical aspects of implementing it on a quantum
device. Specifically, we compute the computational cost of implementing the
two Ansätze in terms of rank, number of variational parameters, the system
size, and we assess the depth and number of variational parameters necessary
to achieve a desired accuracy. Furthermore, we evaluate the algorithms
performance in noisy quantum hardware when simulating the same physical model.

The article is structured as follows. In Section 2, we review the theoretical
background for Markovian open quantum systems and variational quantum time
evolution, and we introduce the quantum algorithm with the two low-rank
parametrized representations of the quantum state. In Section 3, we benchmark
and compare the performance of two different Ansätze, both within statevector
simulation and with Qiskit’s Qasm simulator including noise. Section 4
presents our conclusions and the outlook of our work.

##  2 Model and methods

The dynamics of an open quantum system interacting with a Markovian
environment is governed by the Lindblad master equation [1]

$\displaystyle\dot{\hat{\rho}}$
$\displaystyle=\mathcal{L}\left[\hat{\rho}\right]=-i\left[\hat{H},\hat{\rho}\right]$
$\displaystyle+\sum_{k}\gamma_{k}\left(\hat{c}_{k}\hat{\rho}\hat{c}_{k}^{{\dagger}}-\frac{1}{2}\left\\{\hat{c}_{k}^{{\dagger}}\hat{c}_{k},\hat{\rho}\right\\}\right)\,,$
(1)

where $\mathcal{L}$ is the Liouvillian superoperator, $\hat{\rho}$ is the
density matrix of the system, $\hat{H}$ is the system’s Hamiltonian,
$\hat{c}_{k}$ is the $k-th$ jump operator which describes the mechanism of
dissipation, and $\gamma_{k}$ is the corresponding dissipation rate. For most
systems of physical interest, both $\hat{H}$ and $\hat{c}_{k}$ are quasi-local
operators, i.e., they can be expressed as sums of low-weight Pauli strings,
whose number scales polynomially with the system size $n$ [42].

###  2.1 Low-rank diagonal Ansatz

Figure 1: Encoding the variational Ansatz in a quantum processor. (a) The
computational basis states $|\mathbf{x}_{i}\rangle$ are prepared in a quantum
circuit by application of Pauli $\hat{X}$ gates. (b) Each state
$|\mathbf{x}_{i}\rangle$ is then rotated according to a sequence of
parameterized gates $\hat{U}$, generating the pure states $|\psi_{i}\rangle$
which enter the statistical mixture. (c) The variational Ansatz
$\hat{\rho}(\bm{(}\beta))$ is expressed as a linear combination of pure
states, where each computational basis state $|\mathbf{x}_{i}\rangle$ is
rotated by the same unitary $\hat{U}(\bm{\theta})$ \- Ansatz I, or
individually by $\hat{U}(\bm{\theta}^{(i)})$ \- Ansatz II.

A density matrix can always be expressed in diagonal form as

$\hat{\rho}=\sum_{i=1}^{N}p_{i}\ket{\psi_{i}}\bra{\psi_{i}}\,,$ (2)

where $\ket{\psi_{i}}$ are a set of not necessarily orthogonal but linearly
independent states. In this work, we assume that the open quantum system is
characterized by low entropy along its dynamics, it is convenient to
approximate it in terms of a truncated-rank Ansatz as [36, 43, 44]

$\hat{\rho}(\bm{\beta})=\sum_{i=1}^{R}\alpha_{i}\hat{U}(\bm{\theta})\ket{\mathbf{x}_{i}}\bra{\mathbf{x}_{i}}\hat{U}(\bm{\theta})^{{\dagger}}\,,$
(3)

where $\mathbf{x}_{i}$ are a set of states selected among the $2^{n}$
computational basis states of $n$ qubits,
$\bm{\beta}=(\bm{\alpha},\bm{\theta})$ are variational parameters,
$\alpha_{i}\in\mathbb{R}^{+}$, $\hat{U}(\bm{\theta})$ is a parameterized
unitary operator which is encoded as a variational quantum circuit. This
ansatz does not assume $\mathrm{Tr}(\hat{\rho})=1$, which must be taken into
account when estimating expectation values of operators. For this variational
Ansatz to be efficient, we need the rank $R$ to be $O(\mbox{poly}(n))$. This
condition is often fulfilled in many systems of physical interest [45, 39, 46,
47, 41]. A possible, and arguably more expressive generalization of Ansatz
(3), consists in introducing a different set of values for the parameters
$\bm{\theta}$ for each state in the mixture, resulting in

$\hat{\rho}(\bm{\beta})=\sum_{i=1}^{R}\alpha_{i}\hat{U}(\bm{\theta^{(i)}})\ket{\mathbf{x}_{i}}\bra{\mathbf{x}_{i}}\hat{U}(\bm{\theta^{(i)}})^{{\dagger}}\,.$
(4)

In what follows, we refer to (3) and (4) as Ansatz I and Ansatz II
respectively. Notice that, while in Ansatz I the states of the statistical
mixture are mutually orthogonal, this is not true in Ansatz II. The two
Ansätze are summarized in the diagram presented in Fig. 1. Here we adopt, for
the parameterized circuit $\hat{U}$, a hardware efficient parametrized circuit
[13, 15, 14], with a number of layers $N_{l}$. The expressive power of Ansatz
I or II depends on both $N_{l}$ and $R$, which measure respectively the amount
of entanglement and the mixedness of the quantum state. It has been suggested
that in physically relevant systems there is a trade-off between these two
features, as intuitively the influence of the environment partially suppresses
the entanglement of the system. Here, however, we will study the effectiveness
of Ansatz I and II as a function of the two parameters independently. Ansatz I
and II have different computational costs, with Ansatz II requiring a larger
number of parameters and measurements in the estimate of expectation values.
In Appendix C, we analyze the resource requirements of the present algorithm
and show that the number of quantum circuits to be measured scales as
$O(R^{2}N_{l}+RN^{2}_{l})$ and $O(R^{2}N_{l}^{2})$ for Ansatz I and II,
respectively.

###  2.2 Mclachlan’s variational principle

We review the essential aspects of Mclachlan’s variational principle [48].
Assume a general variational quantum state $\hat{\rho}(\bm{\beta}(t))$, where
$\bm{\beta}$ is a set of real-valued parameters. Evolution over a small time
interval $\delta t$ is described by a variation of $\bm{\beta}$ defined by

$\displaystyle\hat{\rho}(\bm{\beta}(t+\delta t))$
$\displaystyle=\hat{\rho}(\bm{\beta}(t))$
$\displaystyle+\sum_{j}\frac{\partial\hat{\rho}(\bm{\beta})}{\partial\beta_{j}}\dot{\beta_{j}}\delta
t+O(\delta t^{2})\,.$ (5)

The same evolution, according to the master equation Eq. (1), is instead given
by

$\displaystyle\hat{\rho}(\bm{\beta}(t+\delta t))$
$\displaystyle=\hat{\rho}(\bm{\beta}(t))$
$\displaystyle+\mathcal{L}\left[\hat{\rho}(\bm{\beta}(t))\right]\delta
t+O(\delta t^{2})\,.$ (6)

The McLachlan’s variational principle follows from the minimization of the
$L^{2}$-distance between states (5) and (6). To leading order in $\delta t$,
this corresponds to minimizing the norm
$\left\lVert\sum_{j}\frac{\partial\hat{\rho}(\bm{\beta})}{\partial\beta_{j}}\dot{\beta_{j}}-\mathcal{L}\left[\hat{\rho}(\bm{\beta})\right]\right\rVert$
with respect to $\dot{\bm{\beta}}$. This leads to the set of differential
equations

$\displaystyle\sum_{j}M_{kj}\dot{\beta_{j}}=V_{k}\,,$ (7) $\displaystyle
M_{kj}=\mathrm{Tr}\left[\frac{\partial\hat{\rho}(\bm{\beta})}{\partial\beta_{k}}\
\frac{\partial\hat{\rho}(\bm{\beta})}{\partial\beta_{j}}\right]\,,$
$\displaystyle
V_{k}=\mathrm{Tr}\left[\frac{\partial\hat{\rho}(\bm{\beta})}{\partial\beta_{k}}\
\mathcal{L}\left[\hat{\rho}(\bm{\beta})\right]\right],$

where $M_{kj}$ are the elements of a positive semi-definite matrix $M$, and
$V_{k}$ the components of an array $V$. The variational quantum time evolution
algorithm can be summarized as follows. For each time step, the entries of $M$
and $V$ are computed by estimating terms of the form

$\displaystyle\bra{\psi_{i}}\hat{O}\ket{\psi_{j}}\,,\leavevmode\nobreak\
\bra{\partial\psi_{i}/\partial\theta_{k}}\hat{O}\ket{\psi_{j}}\,,$
$\displaystyle\bra{\psi_{i}}\hat{O}\ket{\partial\psi_{j}/\partial\theta_{k}}\,,\leavevmode\nobreak\
\bra{\partial\psi_{i}/\partial\theta_{k}}\hat{O}\ket{\partial\psi_{j}/\partial\theta_{k}}\,,$

on a quantum device, where $\hat{O}$ is a unitary operator typically given by
a low-weight Pauli string, and
$\ket{\psi_{i}}=\hat{U}(\bm{\beta})\ket{\mathbf{x}_{i}}$ and
$\ket{\psi_{i}}=\hat{U}(\bm{\beta}^{(i)})\ket{\mathbf{x}_{i}}$, for Ansatz I
and II, respectively. Then, the linear system Eq. (7) is solved on a classical
processor and the variational parameters are updated. This step is iterated
until the final time is reached. The expressions of Eqs. (7) and of the
quantum circuits are detailed in Sections A.2 and A.3 of the Appendix. Note
that the present formulation is not trace-preserving. Previous studies [49]
have addressed this issue by enforcing trace conservation with Lagrange
multipliers. Here, we adopt an approach consisting in estimating expectation
values of observables $\hat{O}$ as

$\langle\hat{O}\rangle=\frac{\mathrm{Tr}[\hat{\rho}(\bm{\beta})\hat{O}]}{\mathrm{Tr}[\hat{\rho}(\bm{\beta})]}=\frac{\sum_{i}\alpha_{i}\bra{\psi_{i}}\hat{O}\ket{\psi_{i}}}{\sum_{i}\alpha_{i}}\,.$
(8)

Finally, it is important to notice that the present algorithm requires only
$n$ qubits, i.e., as many as are needed to encode the pure states
$\ket{\psi_{i}}$. This comes as a considerable advantage compared to most
algorithms relying on the purification of the quantum state [48] or
vectorization of the density matrix [27], which instead require $2n$ qubits.
This advantage is particularly relevant in view of the application of error
mitigation strategies [50, 51, 52, 53, 54, 55, 56] whose effectiveness depends
critically both on the width and depth of a quantum circuit.

##  3 Numerical Results

In this section, we demonstrate the low-rank quantum time evolution (LRQTE)
algorithm on the dissipative TFIM. This model is characterized by the
Hamiltonian

$\hat{H}_{TFI}=J_{z}\sum_{\langle
j,k\rangle}\hat{\sigma}_{j}^{z}\hat{\sigma}_{k}^{z}+h\sum_{j}\hat{\sigma}_{j}^{x}\,,$
(9)

and by the dissipative dynamics generated by the jump operators
$\hat{c}_{k}=\hat{\sigma}_{k}^{-}=\frac{\hat{\sigma}_{k}^{x}-i\hat{\sigma}_{k}^{y}}{2}$.
Here, $\hat{\sigma}_{j}^{\alpha}$ are Pauli matrices and we assume a lattice
of interacting spin-1/2 sites, where $J_{z}$ is the coupling strength, $h$ is
an external magnetic field, and $\langle j,k\rangle$ indexes pairs of nearest
neighbouring sites. We consider the case with $J_{z}=1$, $h=0.5$, and
$\gamma_{k}=\gamma=1$. We further assume the initial state to be
$\ket{\downarrow\downarrow...\downarrow}\bra{\downarrow\downarrow...\downarrow}$,
i.e., the all-spin down configuration. With this choice of system parameters
and initial state, the Lindblad dynamics towards the steady state is low-rank,
as will appear from the simulations. As a proof-of-concept, we present the
numerical results of the noise-free statevector simulation with Julia’s Yao
package [57] of a $3\times 3$ square lattice with open boundary conditions,
which we use as a benchmark for the algorithm. Specifically, we compute the
average magnetization along the $z$-axis
$s_{z}=\mathrm{Tr}[\hat{\rho}\hat{\sigma}_{z}]$, with
$\hat{\sigma}_{z}=\sum_{i=1}^{n}\hat{\sigma}_{i}^{z}/n$, the average
magnetization along the $x$-axis
$s_{x}=\mathrm{Tr}[\hat{\rho}\hat{\sigma}_{x}]$, with
$\hat{\sigma}_{x}=\sum_{i=1}^{n}\hat{\sigma}_{i}^{x}/n$, the purity of the
variational state $P=\mathrm{Tr}\left[\hat{\rho}^{2}\right]$, and the
infidelity $1-\mathcal{F}$, where
$\mathcal{F}=\mathrm{Tr}\left[\sqrt{\sqrt{\hat{\rho}}\hat{\rho}_{0}\sqrt{\hat{\rho}}}\right]^{2}$
is the fidelity computed with respect to the exact time-evolved state
$\hat{\rho}_{0}$ obtained with the package QuTiP [58, 59] from the direct
integration of the Lindblad master equation. We consider here a hardware-
efficient parametrized quantum circuit $\hat{U}(\bm{\theta})$ inspired by the
Trotter-Suzuki decomposition of $e^{-idt\hat{H}_{TFI}}$ [60]

$\begin{split}\hat{U}(\bm{\theta})&=\prod_{l=1}^{N_{l}}\hat{U}_{l}(\theta_{l})\\\
&=\prod_{l=1}^{N_{l}}\left[\prod_{j}e^{-i\hat{\sigma}_{j}^{x}\frac{\theta_{j,l}}{2}}\right]\left[\prod_{\langle
j,k\rangle}e^{-i\hat{\sigma}_{j}^{z}\hat{\sigma}_{k}^{z}\frac{\theta_{\langle
j,k\rangle,l}}{2}}\right]\,.\end{split}$ (10)

We compare the accuracy and the resource efficiency of the two variational
Ansätze, and investigate the importance of the choice of the truncated basis
for the LRQTE. In addition, a simulation of a smaller instance of the one-
dimensional dissipative TFIM is carried out on Qiskit’s Qasm Simulator [61],
accounting for noise and the connectivity of IBM’s Lagos quantum device.

###  3.1 Statevector Simulation

Figure 2: Time evolution of a $3\times 3$ dissipative TFIM with
$J_{z}=\gamma=1$, $h=0.5$ and $dt=0.005$, simulated with Ansatz I with
$\textrm{rank}=10$: (a) average magnetization along the $x$-axis $s_{x}$, (b)
average magnetization along the $z$-axis $s_{z}$, (c) purity of the
variational state $\hat{\rho}$, (d) infidelity $1-\mathcal{F}$ between the
variational and exact quantum states.

In Figs. 2 and 4 we present the results of the noise-free statevector
simulation using Ansatz I and II, respectively.

Figure 3: Infidelity between the variational and exact quantum states,
computed with Ansatz I with 10 layers and rank $R=\\{10,20,30,46\\}$. Figure
4: Time evolution of a $3\times 3$ dissipative TFIM with $J_{z}=\gamma=1$,
$h=0.5$ and $dt=5\times 10^{-3}$ computed with Ansatz II with rank $R=10$: (a)
average magnetization along the $x$-axis $s_{x}$, (b) average magnetization
along the $z$-axis $s_{z}$, (c) purity of the variational state $\hat{\rho}$,
(d) infidelity $1-\mathcal{F}$ between the variational and exact quantum
states.

Overall, both Ansätze achieve very a high fidelity, highlighting the
effectiveness of the algorithm. In particular, both schemes excel at capturing
the final steady-state. Notably, the infidelity displays a maximum around
$\gamma t=2$, which coincides with a minimum in the state purity. For Ansatz
I, the peak infidelity decreases when increasing the rank $R$ (see Fig. 3),
but not when increasing the number of layers (see Fig. 2(d)). This result
indicates that Ansatz I is less effective in representing faithfully all the
distinct states occurring in the statistical mixture. Conversely, increasing
the number of layers of Ansatz II is sufficient to decrease the overall
infidelity. Finally, we compare the integrated Bures distance [62] to the
exact solution, given by

$I_{B}(T)=\frac{1}{T}\int_{0}^{T}\sqrt{2\textrm{Tr}\left[\sqrt{\hat{\rho}}\hat{\rho}_{0}\sqrt{\hat{\rho}}\right]}dt\,,$
(11)

where, $\hat{\rho}_{0}$ is the exact state of the system, computed with QuTiP,
at $\gamma t=7$. The results are presented in Fig. 5, which illustrates the
trade-off between the computational resources and accuracy. As expected,
Ansatz II achieves the same accuracy as Ansatz I with fewer layers. However,
it requires a considerably larger number of parameters and circuit
measurements to achieve the same target accuracy. Ultimately, the choice
between the two Ansätze depends on the problem at hand. Specifically, whether
we need to optimize time and memory requirements, or rather keep a low circuit
depth to minimize quantum error. In Appendix C, a derivation of the
computational cost of both Ansätze is presented, in terms of the number of
distinct quantum circuits required to estimate the entries of $M$ and $V$.

Figure 5: The accuracy measured by the integrated Bures error $I_{B}$ at
$\gamma t=7$ as a function of (a) the number of layers and (b) the total
number of angular variational parameters for Ansatz I and II.

###  3.2 Simulation with noise

Figure 6: Time evolution of a 2-site dissipative TFIM with $J_{z}=\gamma=1$,
$h=0.5$ and $dt=1\times 10^{-2}$ computed with Ansatz I with rank $R=4$: (a)
average magnetization along the $x$-axis $s_{x}$, (b) average magnetization
along the $z$-axis $s_{z}$, (c) purity of the variational state $\hat{\rho}$,
(d) infidelity $1-\mathcal{F}$ between the variational and exact quantum
states. Each curve was obtained with different quantum simulators: ’Exact’ -
QuTiP, ’SV’ - statevector, ’Ideal’ - IBM’s Qasm Simulator with no hardware
noise, and ’Noisy’ - Qasm Simulator, with the noise model of IBM’s Lagos
machine. For the ’Ideal’ and ’Noisy’ curves, each point in time is an average
of 10 simulations of the corresponding time step and the shaded area is their
standard deviation. For the Qiskit simulations, each quantum circuit is run
20000 times.

We simulate a 2-site and a 3-site (1D) dissipative TFIM model using IBM’s Qasm
Simulator [61]. We restrict this analysis to Ansatz I and $N_{l}=2$, which
represents faithfully the state of such small systems. Results are displayed
in Fig. 6. We compare results obtained using a noise model corresponding to
the IBM Lagos quantum device, with noiseless simulations (i.e., including only
the shot noise from projective readout), and with statevector simulations.
Noisy estimates of the matrix elements can lead to a non semi-positive defined
matrix $M$. As the matrix is typically singular and must be regularized to
solve the system Eq. (7), the occurrence of negative eigenvalues calls for a
different regularization scheme, as detailed in Appendix F. For details on the
estimation of $M$ and $V$, we refer the reader to Appendix D. Fig.7 presents
the results of the full-rank QTE for the 2-site dissipative TFIM. We compare
the results of the statevector, noiseless, and noisy simulations of the $x$\-
and $z$-polarization, and we study the purity and the fidelity to the exact
results. To compute the purity and the fidelity in the noiseless and noisy
simulations, a full estimate of the density matrix at each time would be
needed. Instead of simulating a computationally demanding full-state
tomography, we resort to a hybrid approach consisting in computing the
variational parameters from the noiseless and noisy simulations, and using
them to compute the density matrix directly from Eq. (3). The results for the
average polarizations indicate that the noiseless simulation reproduces
accurately the statevector simulation, when using the same regularization
($\lambda=0.04$, see Appendix F for details about the regularization parameter
$\lambda$). However, the noisy simulation shows a significant bias, which is
due to both the limited accuracy in estimating $M$ and $V$, and the error in
the estimate of the observables. This analysis is further supported by the
study of the purity and infidelity. We also demonstrate the algorithm on a
3-site dissipative TFIM with open boundary conditions, taking $R=3$ and
selecting the computational basis vectors
$\\{\ket{000},\ket{001},\ket{010}\\}$ (Fig. 7). The statevector simulation
deviates from the exact result due to the low-rank approximation, as confirmed
by the statevector simulation with $R=8$. As in the previous case, the noisy
simulation results in a significant bias for the observables. The low-rank
approximation also affects the purity and fidelity. In particular, now the
infidelity stays finite in the long-time limit.

Figure 7: Time evolution of a 3-site dissipative TFIM $J_{z}=\gamma=1$,
$h=0.5$ and $dt=1\times 10^{-2}$ computed with Ansatz I with rank $R=3$ (solid
lines) and $8$ (dash-dotted lines): (a) average magnetization along the
$x$-axis $s_{x}$, (b) average magnetization along the $z$-axis $s_{z}$, (c)
purity of the variational state $\hat{\rho}$, (d) infidelity $1-\mathcal{F}$
between the variational and exact quantum states. For the ’Ideal’ and ’Noisy’
curves, each point in time is an average of 10 simulations of the
corresponding time step and the shaded area is the standard deviation. For the
Qiskit simulations, each quantum circuit is run 20000 times.

##  4 Conclusions & Outlook

In this paper, we propose a novel time-dependent variational quantum
algorithm, which directly integrates the Lindblad master equation. In this
framework, the state of the system is approximated by a low-rank Ansatz,
defined as a linear combination of pure states that are represented by
parameterized quantum circuits. The time-evolution of the system is then
achieved by evolving the variational parameters, according to Mclachlan’s
variational principle. This algorithm requires only evaluating circuits with
$n$ qubits, less than other existing similar algorithms.

As a demonstration of the algorithm and a proof of concept, we applied it to a
dissipative transverse field Ising system. Two types of simulations were
performed, one based on the measurement of quantum circuits with Qiskit’s Qasm
Simulator, with and without hardware noise, and a noise-free statevector
simulation, which allowed us to verify the performance of the algorithm when
simulating larger systems and study the expressivity error of the variational
Ansätze.

The statevector simulations were implemented with two different Ansätze,
respectively expressing a statistical mixture of orthogonal or non-orthogonal
pure states. We find that the orthogonal Ansatz is more resource-efficient.
Nevertheless, the non-orthogonal scheme has a more flexible structure and can
achieve better accuracy with shallower quantum circuits.

Statevector simulations showed that the algorithm can efficiently model the
open system dynamics in cases where the system entropy remains bounded, at a
computational cost which scales only linearly with the rank.

A considerable advantage of the present algorithm is that it requires only $n$
qubits, i.e., as many as are needed to encode the pure states
$\ket{\psi_{i}}$. This comes as a considerable advantage compared to most
algorithms relying on purification of the quantum state [48] or vectorization
of the density matrix [27], which instead require $2n$ qubits. This advantage
is particularly relevant in view of the application of error mitigation
strategies [50, 51, 52, 53, 54, 55, 56] whose effectiveness depends critically
both on the width and depth of a quantum circuit.

Simulations accounting for noise on Qiskit’s Qasm simulator showed that
accurate regularization of the matrix $M$ defining the variational dynamics is
needed in order to avoid numerical instabilities. Here, we applied the
diagonal shift regularization method, which preserves the positivity of $M$.
Noise may additionally lead to biases in the estimate of expectation values,
hence to inaccurate evolution of the variational parameters. This issue can be
at least partially addressed by quantum error mitigation techniques, which
will demand further investigation.

A possible generalization of the present algorithm is the implementation of an
adaptive-rank scheme, similarly to what has been proposed for classical
simulations of open quantum dynamics [41, 39, 46, 47, 45]. An adaptive-rank
scheme would avoid accounting for irrelevant states, i.e., with vanishing
probabilities, in the statistical mixture, thereby increasing computational
efficiency while reducing numerical instabilities due to small probabilities.
An adaptive-rank scheme may also be used in cases where the transient dynamics
is not necessarily low-rank while the long-time limit is. In these cases, an
adaptive bounded-rank simulation of a fictitious system dynamics would
efficiently describe the long-time limit.

Another possible improvement is the generalization of the projected
variational quantum dynamics to the present low-rank algorithm. For unitary
quantum dynamics, the p-VQD approach was successful in addressing the
occurrence of biases and uncontrolled variance in Monte-Carlo estimates,
arising from vanishing amplitudes in the variational wave functions [60, 63,
64]. We expect a similar improvement in the present case.

##  5 Acknowledgements

Sara Santos and Xinyu Song contributed equally to this work. We thank Zakari
Denis for insightful discussions. We acknowledge the use of IBM Quantum
services for this work. The views expressed are those of the authors, and do
not reflect the official policy or position of IBM or the IBM Quantum team.
This work was supported by the Swiss National Science Foundation through
Projects No. 200020_185015, 200020_215172, and 20QU-1_215928, and was
conducted with the financial support of the EPFL Science Seed Fund 2021.

## References

  * Breuer and Petruccione [2007] Heinz-Peter Breuer and Francesco Petruccione.  _The Theory of Open Quantum Systems_.  Oxford University Press, 2007.
  * Daley [2014] Andrew J. Daley.  Quantum trajectories and open many-body quantum systems.  _Advances in Physics_ , 63(2):77–149, mar 2014.  doi: 10.1080/00018732.2014.933502.  URL https://doi.org/10.1080%2F00018732.2014.933502.
  * Dum et al. [1992] R. Dum, P. Zoller, and H. Ritsch.  Monte carlo simulation of the atomic master equation for spontaneous emission.  _Phys. Rev. A_ , 45:4879–4887, Apr 1992.  doi: 10.1103/PhysRevA.45.4879.  URL https://link.aps.org/doi/10.1103/PhysRevA.45.4879.
  * Dalibard et al. [1992] Jean Dalibard, Yvan Castin, and Klaus Mølmer.  Wave-function approach to dissipative processes in quantum optics.  _Phys. Rev. Lett._ , 68:580–583, Feb 1992.  doi: 10.1103/PhysRevLett.68.580.  URL https://link.aps.org/doi/10.1103/PhysRevLett.68.580.
  * Charmichael [1993] Howard Charmichael.  _An Open Systems Approach to Quantum Optics_.  Springer-Verlag, 1993.  ISBN 978-0-387-56634-4.
  * Weimer et al. [2021] Hendrik Weimer, Augustine Kshetrimayum, and Román Orús.  Simulation methods for open quantum many-body systems.  _Rev. Mod. Phys._ , 93:015008, Mar 2021.  doi: 10.1103/RevModPhys.93.015008.  URL https://link.aps.org/doi/10.1103/RevModPhys.93.015008.
  * Preisser et al. [2023] Guillermo Preisser, David Wellnitz, Thomas Botzung, and Johannes Schachenmayer.  Comparing bipartite entropy growth in open-system matrix-product simulation methods.  _Phys. Rev. A_ , 108:012616, Jul 2023.  doi: 10.1103/PhysRevA.108.012616.  URL https://link.aps.org/doi/10.1103/PhysRevA.108.012616.
  * Feynman [1982] Richard P. Feynman.  _Simulating physics with computers_.  International Journal of Theoretical Physics, 1982.  doi: 10.1007/BF02650179.
  * Lloyd [1996] Seth Lloyd.  Universal quantum simulators.  _Science_ , 273(5278):1073–1078, 1996.  doi: 10.1126/science.273.5278.1073.
  * Tacchino et al. [2020] Francesco Tacchino, Alessandro Chiesa, Stefano Carretta, and Dario Gerace.  Quantum computers as universal quantum simulators: State-of-the-art and perspectives.  _Advanced Quantum Technologies_ , 3(3):1900052, 2020.  doi: https://doi.org/10.1002/qute.201900052.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/qute.201900052.
  * Childs et al. [2021] Andrew M. Childs, Yuan Su, Minh C. Tran, Nathan Wiebe, and Shuchen Zhu.  Theory of trotter error with commutator scaling.  _Phys. Rev. X_ , 11:011020, Feb 2021.  doi: 10.1103/PhysRevX.11.011020.  URL https://link.aps.org/doi/10.1103/PhysRevX.11.011020.
  * Preskill [2018] John Preskill.  Quantum Computing in the NISQ era and beyond.  _Quantum_ , 2:79, August 2018.  ISSN 2521-327X.  doi: 10.22331/q-2018-08-06-79.  URL https://doi.org/10.22331/q-2018-08-06-79.
  * Cerezo et al. [2021] M. Cerezo, A. Arrasmith, R. Babbush, and et al.  Variational quantum algorithms.  _Nat Rev Phys_ , 3:625–644, 2021.  doi: https://doi.org/10.1038/s42254-021-00348-9.
  * Bauer et al. [2020] Bela Bauer, Sergey Bravyi, Mario Motta, and Garnet Kin-Lic Chan.  Quantum algorithms for quantum chemistry and quantum materials science.  _Chemical Reviews_ , 120(22):12685–12717, 2020.  doi: 10.1021/acs.chemrev.9b00829.  URL https://doi.org/10.1021/acs.chemrev.9b00829.  PMID: 33090772.
  * Bharti et al. [2022] Kishor Bharti, Alba Cervera-Lierta, Thi Ha Kyaw, Tobias Haug, Sumner Alperin-Lea, Abhinav Anand, Matthias Degroote, Hermanni Heimonen, Jakob S. Kottmann, Tim Menke, Wai-Keong Mok, Sukin Sim, Leong-Chuan Kwek, and Alán Aspuru-Guzik.  Noisy intermediate-scale quantum algorithms.  _Rev. Mod. Phys._ , 94:015004, Feb 2022.  doi: 10.1103/RevModPhys.94.015004.  URL https://link.aps.org/doi/10.1103/RevModPhys.94.015004.
  * Barreiro et al. [2011] Julio T. Barreiro, Markus Müller, Philipp Schindler, Daniel Nigg, Thomas Monz, Michael Chwalla, Markus Hennrich, Christian F. Roos, Peter Zoller, and Rainer Blatt.  An open-system quantum simulator with trapped ions.  _Nature_ , 470(7335):486–491, February 2011.  ISSN 1476-4687.  doi: 10.1038/nature09801.  URL http://dx.doi.org/10.1038/nature09801.
  * Wang et al. [2011] Hefeng Wang, S. Ashhab, and Franco Nori.  Quantum algorithm for simulating the dynamics of an open quantum system.  _Phys. Rev. A_ , 83:062317, Jun 2011.  doi: 10.1103/PhysRevA.83.062317.  URL https://link.aps.org/doi/10.1103/PhysRevA.83.062317.
  * Sweke et al. [2014] Ryan Sweke, Ilya Sinayskiy, and Francesco Petruccione.  Simulation of single-qubit open quantum systems.  _Phys. Rev. A_ , 90:022331, Aug 2014.  doi: 10.1103/PhysRevA.90.022331.  URL https://link.aps.org/doi/10.1103/PhysRevA.90.022331.
  * Sweke et al. [2015] Ryan Sweke, Ilya Sinayskiy, Denis Bernard, and Francesco Petruccione.  Universal simulation of markovian open quantum systems.  _Phys. Rev. A_ , 91:062308, Jun 2015.  doi: 10.1103/PhysRevA.91.062308.  URL https://link.aps.org/doi/10.1103/PhysRevA.91.062308.
  * Di Candia et al. [2015] R. Di Candia, J. S. Pedernales, A. del Campo, E. Solano, and J. Casanova.  Quantum simulation of dissipative processes without reservoir engineering.  _Scientific Reports_ , 5(1), May 2015.  ISSN 2045-2322.  doi: 10.1038/srep09981.  URL http://dx.doi.org/10.1038/srep09981.
  * Sweke et al. [2016] R. Sweke, M. Sanz, I. Sinayskiy, F. Petruccione, and E. Solano.  Digital quantum simulation of many-body non-markovian dynamics.  _Phys. Rev. A_ , 94:022317, Aug 2016.  doi: 10.1103/PhysRevA.94.022317.  URL https://link.aps.org/doi/10.1103/PhysRevA.94.022317.
  * Chenu et al. [2017] A. Chenu, M. Beau, J. Cao, and A. del Campo.  Quantum simulation of generic many-body open system dynamics using classical noise.  _Phys. Rev. Lett._ , 118:140403, Apr 2017.  doi: 10.1103/PhysRevLett.118.140403.  URL https://link.aps.org/doi/10.1103/PhysRevLett.118.140403.
  * Cleve and Wang [2017] Richard Cleve and Chunhao Wang.  Efficient Quantum Algorithms for Simulating Lindblad Evolution.  In Ioannis Chatzigiannakis, Piotr Indyk, Fabian Kuhn, and Anca Muscholl, editors, _44th International Colloquium on Automata, Languages, and Programming (ICALP 2017)_ , volume 80 of _Leibniz International Proceedings in Informatics (LIPIcs)_ , pages 17:1–17:14, Dagstuhl, Germany, 2017. Schloss Dagstuhl – Leibniz-Zentrum für Informatik.  ISBN 978-3-95977-041-5.  doi: 10.4230/LIPIcs.ICALP.2017.17.  URL https://drops.dagstuhl.de/entities/document/10.4230/LIPIcs.ICALP.2017.17.
  * García-Pérez et al. [2020] Guillermo García-Pérez, Matteo A. C. Rossi, and Sabrina Maniscalco.  Ibm q experience as a versatile experimental testbed for simulating open quantum systems.  _npj Quantum Information_ , 6(1), January 2020.  ISSN 2056-6387.  doi: 10.1038/s41534-019-0235-y.  URL http://dx.doi.org/10.1038/s41534-019-0235-y.
  * Su and Li [2020a] Hong-Yi Su and Ying Li.  Quantum algorithm for the simulation of open-system dynamics and thermalization.  _Phys. Rev. A_ , 101:012328, Jan 2020a.  doi: 10.1103/PhysRevA.101.012328.  URL https://link.aps.org/doi/10.1103/PhysRevA.101.012328.
  * Endo et al. [2020] Suguru Endo, Jinzhao Sun, Ying Li, Simon C. Benjamin, and Xiao Yuan.  Variational quantum simulation of general processes.  _Phys. Rev. Lett._ , 125:010501, Jun 2020.  doi: 10.1103/PhysRevLett.125.010501.  URL https://link.aps.org/doi/10.1103/PhysRevLett.125.010501.
  * Yoshioka et al. [2020] Nobuyuki Yoshioka, Yuya O. Nakagawa, Kosuke Mitarai, and Keisuke Fujii.  Variational quantum algorithm for nonequilibrium steady states.  _Phys. Rev. Res._ , 2:043289, November 2020.  doi: 10.1103/PhysRevResearch.2.043289.  URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.043289.
  * Ramusat and Savona [2021] Nathan Ramusat and Vincenzo Savona.  A quantum algorithm for the direct estimation of the steady state of open quantum systems.  _Quantum_ , 5:399, February 2021.  ISSN 2521-327X.  doi: 10.22331/q-2021-02-22-399.  URL https://doi.org/10.22331/q-2021-02-22-399.
  * Han et al. [2021] J. Han, W. Cai, L. Hu, X. Mu, Y. Ma, Y. Xu, W. Wang, H. Wang, Y. P. Song, C.-L. Zou, and L. Sun.  Experimental simulation of open quantum system dynamics via trotterization.  _Phys. Rev. Lett._ , 127:020504, Jul 2021.  doi: 10.1103/PhysRevLett.127.020504.  URL https://link.aps.org/doi/10.1103/PhysRevLett.127.020504.
  * Kamakari et al. [2022] Hirsh Kamakari, Shi-Ning Sun, Mario Motta, and Austin J. Minnich.  Digital quantum simulation of open quantum systems using quantum imaginary–time evolution.  _PRX Quantum_ , 3:010320, February 2022.  doi: 10.1103/PRXQuantum.3.010320.  URL https://link.aps.org/doi/10.1103/PRXQuantum.3.010320.
  * Liu et al. [2022] Zidu Liu, L.-M. Duan, and Dong-Ling Deng.  Solving quantum master equations with deep quantum neural networks.  _Phys. Rev. Res._ , 4:013097, Feb 2022.  doi: 10.1103/PhysRevResearch.4.013097.  URL https://link.aps.org/doi/10.1103/PhysRevResearch.4.013097.
  * Hu et al. [2022] Zixuan Hu, Kade Head-Marsden, David A. Mazziotti, Prineha Narang, and Sabre Kais.  A general quantum algorithm for open quantum dynamics demonstrated with the Fenna-Matthews-Olson complex.  _Quantum_ , 6:726, May 2022.  ISSN 2521-327X.  doi: 10.22331/q-2022-05-30-726.  URL https://doi.org/10.22331/q-2022-05-30-726.
  * Watad and Lindner [2023] Tasneem Watad and Netanel H. Lindner.  Variational quantum algorithms for simulation of lindblad dynamics, May 2023.  URL https://arxiv.org/abs/2305.02815.
  * Chen et al. [2024] Huo Chen, Niladri Gomes, Siyuan Niu, and Wibe Albert de Jong.  Adaptive variational simulation for open quantum systems.  _Quantum_ , 8:1252, February 2024.  ISSN 2521-327X.  doi: 10.22331/q-2024-02-13-1252.  URL https://doi.org/10.22331/q-2024-02-13-1252.
  * Terhal and DiVincenzo [2000] Barbara M. Terhal and David P. DiVincenzo.  Problem of equilibration and the computation of correlation functions on a quantum computer.  _Phys. Rev. A_ , 61:022301, Jan 2000.  doi: 10.1103/PhysRevA.61.022301.  URL https://link.aps.org/doi/10.1103/PhysRevA.61.022301.
  * Ezzell et al. [2023] Nic Ezzell, Elliott M Ball, Aliza U Siddiqui, Mark M Wilde, Andrew T Sornborger, Patrick J Coles, and Zoë Holmes.  Quantum mixed state compiling.  _Quantum Science and Technology_ , 8(3):035001, apr 2023.  doi: 10.1088/2058-9565/acc4e3.  URL https://dx.doi.org/10.1088/2058-9565/acc4e3.
  * Su and Li [2020b] Hong-Yi Su and Ying Li.  Quantum algorithm for the simulation of open-system dynamics and thermalization.  _Phys. Rev. A_ , 101:012328, Jan 2020b.  doi: 10.1103/PhysRevA.101.012328.  URL https://link.aps.org/doi/10.1103/PhysRevA.101.012328.
  * Zhang et al. [2022] Yuxuan Zhang, Shahin Jahanbani, Daoheng Niu, Reza Haghshenas, and Andrew C. Potter.  Qubit-efficient simulation of thermal states with quantum tensor networks.  _Phys. Rev. B_ , 106:165126, Oct 2022.  doi: 10.1103/PhysRevB.106.165126.  URL https://link.aps.org/doi/10.1103/PhysRevB.106.165126.
  * Donatella et al. [2021] Kaelan Donatella, Zakari Denis, Alexandre Le Boité, and Cristiano Ciuti.  Continuous-time dynamics and error scaling of noisy highly entangling quantum circuits.  _Physical Review A_ , 104(6), December 2021.  ISSN 2469-9934.  doi: 10.1103/physreva.104.062407.  URL http://dx.doi.org/10.1103/PhysRevA.104.062407.
  * Cîrstoiu et al. [2020] C. Cîrstoiu, Z. Holmes, J. Iosue, and et al.  Variational fast forwarding for quantum simulation beyond the coherence time.  _npj Quantum Inf_ , 6:82, 2020.  doi: https://doi.org/10.1038/s41534-020-00302-0.
  * Gravina and Savona [2023] Luca Gravina and Vincenzo Savona.  Adaptive variational low-rank dynamics for open quantum systems, November 2023.  URL https://arxiv.org/abs/2312.13676.
  * Kempe et al. [2006] Julia Kempe, Alexei Kitaev, and Oded Regev.  The complexity of the local hamiltonian problem.  _SIAM Journal on Computing_ , 35(5):1070–1097, 2006.  doi: 10.1137/S0097539704445226.  URL https://doi.org/10.1137/S0097539704445226.
  * Verdon et al. [2019] Guillaume Verdon, Jacob Marks, Sasha Nanda, Stefan Leichenauer, and Jack Hidary.  Quantum hamiltonian-based models and the variational quantum thermalizer algorithm, October 2019.  URL https://arxiv.org/abs/1910.02071.
  * Liu et al. [2021] Jin-Guo Liu, Liang Mao, Pan Zhang, and Lei Wang.  Solving quantum statistical mechanics with variational autoregressive networks and quantum circuits.  _Machine Learning: Science and Technology_ , 2(2), February 2021.  doi: 10.1088/2632-2153/aba19d.  URL https://dx.doi.org/10.1088/2632-2153/aba19d.
  * McCaul et al. [2021] Gerard McCaul, Kurt Jacobs, and Denys I. Bondar.  Fast computation of dissipative quantum systems with ensemble rank truncation.  _Physical Review Research_ , 3(1), January 2021.  ISSN 2643-1564.  doi: 10.1103/physrevresearch.3.013017.  URL http://dx.doi.org/10.1103/PhysRevResearch.3.013017.
  * Chen et al. [2021] Yi-Ting Chen, Collin Farquhar, and Robert M. Parrish.  Low-rank density-matrix evolution for noisy quantum circuits.  _npj Quantum Information_ , 7(1):61, Apr 2021.  ISSN 2056-6387.  doi: 10.1038/s41534-021-00392-4.  URL https://doi.org/10.1038/s41534-021-00392-4.
  * Le Bris and Rouchon [2013] C. Le Bris and P. Rouchon.  Low-rank numerical approximations for high-dimensional lindblad equations.  _Physical Review A_ , 87(2), February 2013.  ISSN 1094-1622.  doi: 10.1103/physreva.87.022125.  URL http://dx.doi.org/10.1103/PhysRevA.87.022125.
  * Yuan et al. [2019] Xiao Yuan, Suguru Endo, Qi Zhao, Ying Li, and Simon C. Benjamin.  Theory of variational quantum simulation.  _Quantum_ , 3:191, October 2019.  ISSN 2521-327X.  doi: 10.22331/q-2019-10-07-191.  URL https://doi.org/10.22331/q-2019-10-07-191.
  * Joubert-Doriol and Izmaylov [2015] Loïc Joubert-Doriol and Artur F. Izmaylov.  Problem-free time-dependent variational principle for open quantum systems.  _The Journal of Chemical Physics_ , 142(13), 04 2015.  ISSN 0021-9606.  doi: 10.1063/1.4916384.  URL https://doi.org/10.1063/1.4916384.  134107\.
  * Cai et al. [2023] Zhenyu Cai, Ryan Babbush, Simon C. Benjamin, Suguru Endo, William J. Huggins, Ying Li, Jarrod R. McClean, and Thomas E. O’Brien.  Quantum error mitigation.  _Rev. Mod. Phys._ , 95:045005, Dec 2023.  doi: 10.1103/RevModPhys.95.045005.  URL https://link.aps.org/doi/10.1103/RevModPhys.95.045005.
  * Mari et al. [2021] Andrea Mari, Nathan Shammah, and William J. Zeng.  Extending quantum probabilistic error cancellation by noise scaling.  _Physical Review A_ , 104(5), nov 2021.  doi: 10.1103/physreva.104.052607.  URL https://doi.org/10.1103%2Fphysreva.104.052607.
  * van den Berg et al. [2023] Ewout van den Berg, Zlatko K. Minev, Abhinav Kandala, and Kristan Temme.  Probabilistic error cancellation with sparse pauli–lindblad models on noisy quantum processors.  _Nature Physics_ , 19(8):1116–1121, Aug 2023.  ISSN 1745-2481.  doi: 10.1038/s41567-023-02042-2.  URL https://doi.org/10.1038/s41567-023-02042-2.
  * Endo et al. [2018] Suguru Endo, Simon C. Benjamin, and Ying Li.  Practical quantum error mitigation for near-future applications.  _Physical Review X_ , 8(3), jul 2018.  doi: 10.1103/physrevx.8.031027.  URL https://doi.org/10.1103%2Fphysrevx.8.031027.
  * Endo et al. [2021] Suguru Endo, Zhenyu Cai, Simon C. Benjamin, and Xiao Yuan.  Hybrid quantum-classical algorithms and quantum error mitigation.  _Journal of the Physical Society of Japan_ , 90(3):032001, mar 2021.  doi: 10.7566/jpsj.90.032001.  URL https://doi.org/10.7566%2Fjpsj.90.032001.
  * Nation et al. [2021] Paul D. Nation, Hwajung Kang, Neereja Sundaresan, and Jay M. Gambetta.  Scalable mitigation of measurement errors on quantum computers.  _PRX Quantum_ , 2(4), November 2021.  ISSN 2691-3399.  doi: 10.1103/prxquantum.2.040326.  URL http://dx.doi.org/10.1103/PRXQuantum.2.040326.
  * Yang et al. [2022] Bo Yang, Rudy Raymond, and Shumpei Uno.  Efficient quantum readout-error mitigation for sparse measurement outcomes of near-term quantum devices.  _Physical Review A_ , 106(1), July 2022.  ISSN 2469-9934.  doi: 10.1103/physreva.106.012423.  URL http://dx.doi.org/10.1103/PhysRevA.106.012423.
  * Luo et al. [2020] Xiu-Zhe Luo, Jin-Guo Liu, Pan Zhang, and Lei Wang.  Yao.jl: Extensible, Efficient Framework for Quantum Algorithm Design.  _Quantum_ , 4:341, October 2020.  ISSN 2521-327X.  doi: 10.22331/q-2020-10-11-341.  URL https://doi.org/10.22331/q-2020-10-11-341.
  * Johansson et al. [2012] J.R. Johansson, P.D. Nation, and Franco Nori.  Qutip: An open-source python framework for the dynamics of open quantum systems.  _Computer Physics Communications_ , 183(8):1760–1772, August 2012.  ISSN 0010-4655.  doi: 10.1016/j.cpc.2012.02.021.  URL http://dx.doi.org/10.1016/j.cpc.2012.02.021.
  * Johansson et al. [2013] J.R. Johansson, P.D. Nation, and Franco Nori.  Qutip 2: A python framework for the dynamics of open quantum systems.  _Computer Physics Communications_ , 184(4):1234–1240, April 2013.  ISSN 0010-4655.  doi: 10.1016/j.cpc.2012.11.019.  URL http://dx.doi.org/10.1016/j.cpc.2012.11.019.
  * Barison et al. [2021] Stefano Barison, Filippo Vicentini, and Giuseppe Carleo.  An efficient quantum algorithm for the time evolution of parameterized circuits.  _Quantum_ , 5:512, July 2021.  ISSN 2521-327X.  doi: 10.22331/q-2021-07-28-512.  URL https://doi.org/10.22331/q-2021-07-28-512.
  * et al [2019] Gadi Aleksandrowicz et al.  Qiskit: An Open-source Framework for Quantum Computing.  2019\.  doi: 10.5281/zenodo.2562111.  URL https://doi.org/10.5281/zenodo.2562111.
  * Laha et al. [2021] Aritra Laha, Agrim Aggarwal, and Santosh Kumar.  Random density matrices: Analytical results for mean root fidelity and the mean-square bures distance.  _Physical Review A_ , 104(2), aug 2021.  doi: 10.1103/physreva.104.022438.  URL https://doi.org/10.1103%2Fphysreva.104.022438.
  * Gacon et al. [2024] Julien Gacon, Jannes Nys, Riccardo Rossi, Stefan Woerner, and Giuseppe Carleo.  Variational quantum time evolution without the quantum geometric tensor.  _Phys. Rev. Res._ , 6:013143, Feb 2024.  doi: 10.1103/PhysRevResearch.6.013143.  URL https://link.aps.org/doi/10.1103/PhysRevResearch.6.013143.
  * Sinibaldi et al. [2023] Alessandro Sinibaldi, Clemens Giuliani, Giuseppe Carleo, and Filippo Vicentini.  Unbiasing time-dependent variational monte carlo by projected quantum evolution.  _Quantum_ , 7:1131, October 2023.  ISSN 2521-327X.  doi: 10.22331/q-2023-10-10-1131.  URL http://dx.doi.org/10.22331/q-2023-10-10-1131.
  * Kimura [2017] Tetsuji Kimura.  Explicit description of the zassenhaus formula.  _Progress of Theoretical and Experimental Physics_ , 2017(4), April 2017.  ISSN 2050-3911.  doi: 10.1093/ptep/ptx044.  URL http://dx.doi.org/10.1093/ptep/ptx044.
  * Mitarai and Fujii [2019] Kosuke Mitarai and Keisuke Fujii.  Methodology for replacing indirect measurements with direct measurements.  _Phys. Rev. Res._ , 1:013006, August 2019.  doi: 10.1103/PhysRevResearch.1.013006.  URL https://link.aps.org/doi/10.1103/PhysRevResearch.1.013006.
  * Zoufal et al. [2023] Christa Zoufal, David Sutter, and Stefan Woerner.  Error bounds for variational quantum time evolution.  _Phys. Rev. Appl._ , 20:044059, Oct 2023.  doi: 10.1103/PhysRevApplied.20.044059.  URL https://link.aps.org/doi/10.1103/PhysRevApplied.20.044059.
  * Medvidović and Sels [2023] Matija Medvidović and Dries Sels.  Variational quantum dynamics of two-dimensional rotor models.  _PRX Quantum_ , 4(4), October 2023.  ISSN 2691-3399.  doi: 10.1103/prxquantum.4.040302.  URL http://dx.doi.org/10.1103/PRXQuantum.4.040302.

##  Appendix A Mclachlan’s variational principle

###  A.1 Derivation in the density matrix formalism

The starting point of our derivation is the variational density matrix
$\hat{\rho}(\bm{\beta})$ and the Lindblad master equation
$\dot{\hat{\rho}}=\mathcal{L}[\hat{\rho}]$. Following the Mclachlan’s
variational principle [48], the time evolution of the variational parameters
is obtained from minimizing the $L^{2}$-distance between the variational state
$\hat{\rho}(\bm{\beta}(t))+\sum_{j}\frac{\partial\hat{\rho}(\bm{\beta})}{\partial\beta_{j}}\dot{\beta_{j}}\delta
t$ and the exact time evolved state
$\hat{\rho}(\bm{\beta}(t))+\mathcal{L}\left[\hat{\rho}(\bm{\beta}(t))\right]\delta
t$. This variational principle can be expressed as follows:

$\bm{\dot{\beta}}=\textrm{argmin}\left(\lVert\sum_{j}\frac{\partial\hat{\rho}(\bm{\beta})}{\partial\beta_{j}}\dot{\beta_{j}}-\mathcal{L}\left[\hat{\rho}(\bm{\beta}(t))\right]\rVert_{2}\right)\,,$
(12)

or equivalently,

$\bm{\dot{\beta}}=\textrm{argmin}\left(C(\bm{\dot{\beta}})\right)\,,$ (13)
${}C(\bm{\dot{\beta}})=\mathrm{Tr}\left[\left(\sum_{j}\frac{\partial\hat{\rho}(\bm{\beta})}{\partial\beta_{j}}\dot{\beta_{j}}-\mathcal{L}\left[\hat{\rho}(\bm{\beta}(t))\right]\right)^{\dagger}\left(\sum_{j}\frac{\partial\hat{\rho}(\bm{\beta})}{\partial\beta_{j}}\dot{\beta_{j}}-\mathcal{L}\left[\hat{\rho}(\bm{\beta}(t))\right]\right)\right]\,,$
(14)

where $C(\bm{\dot{\beta}})$ can be interpreted as a cost function. To find the
minimum of $C(\bm{\dot{\beta}})$, we calculate the derivative of
$C(\bm{\dot{\beta}})$ with respect to the variable $\dot{\beta_{k}}$

$\displaystyle\begin{split}\frac{\partial
C(\bm{\dot{\beta}})}{\partial\dot{\beta_{k}}}&=\mathrm{Tr}\left[\frac{\partial\hat{\rho}^{\dagger}}{\partial\beta_{k}}\left(\sum_{j}\frac{\partial\hat{\rho}}{\partial\beta_{j}}\dot{\beta_{j}}-\mathcal{L}[\hat{\rho}]\right)+\left(\sum_{j}\frac{\partial\hat{\rho}}{\partial\beta_{j}}\dot{\beta_{j}}-\mathcal{L}[\hat{\rho}]\right)^{\dagger}\frac{\partial\hat{\rho}}{\partial\beta_{k}}\right]\\\
&=\sum_{j}\left(\mathrm{Tr}\left[\frac{\partial\hat{\rho}^{\dagger}}{\partial{\beta_{k}}}\frac{\partial\hat{\rho}}{\partial\beta_{j}}\right]\dot{\beta_{j}}+h.c.\right)-\left(\mathrm{Tr}\left[\frac{\partial\hat{\rho}^{\dagger}}{\partial\beta_{k}}\mathcal{L}[\hat{\rho}]\right]+h.c.\right)\\\
&=2\left(\sum_{j}\mathcal{R}e\left[\mathrm{Tr}\left[\frac{\partial\hat{\rho}}{\partial{\beta_{k}}}\frac{\partial\hat{\rho}}{\partial\beta_{j}}\right]\right]\dot{\beta_{j}}-\mathcal{R}e\left[\mathrm{Tr}\left[\frac{\partial\hat{\rho}}{\partial{\beta_{k}}}\mathcal{L}[\hat{\rho}]\right]\right]\right)\\\
&=2\left(\sum_{j}\mathrm{Tr}\left[\frac{\partial\hat{\rho}}{\partial{\beta_{k}}}\frac{\partial\hat{\rho}}{\partial\beta_{j}}\right]\dot{\beta_{j}}-\mathrm{Tr}\left[\frac{\partial\hat{\rho}}{\partial{\beta_{k}}}\mathcal{L}[\hat{\rho}]\right]\right)\,,\end{split}$
(15)

where we have used the hermiticity of the density matrix. By setting
$\frac{\partial C(\bm{\dot{\beta}})}{\partial\dot{\beta_{k}}}=0$, we obtain
the equations of motion (EOM) of the variational parameters, presented in Eq.
(7).

###  A.2 Equations of motion - Ansatz I

In this section we derive the EOM of the parameters
$\bm{\beta}=(\bm{\alpha},\bm{\theta})$ by computing the elements $M_{kj}$ and
$V_{k}$ according to Eq. (7). Let
$\hat{\rho}(\bm{\alpha},\bm{\theta})=\sum_{p=1}^{R}\alpha_{p}\hat{U}(\bm{\theta})\ket{\mathbf{x}_{p}}\bra{\mathbf{x}_{p}}\hat{U}(\bm{\theta})^{\dagger}$
be the density matrix and $\ket{\mathbf{x}_{p}}$ computational basis states.
For simplicity, we use the notation $\hat{U}(\bm{\theta})=\hat{U}$. We begin
by writing the derivatives:

$\frac{\partial\hat{\rho}}{\partial\alpha_{p}}=\hat{U}\ket{\mathbf{x}_{p}}\bra{\mathbf{x}_{p}}\hat{U}^{\dagger}\,,$
(16)

and

$\frac{\partial\hat{\rho}}{\partial\theta_{j}}=\sum_{p}\alpha_{p}\left(\frac{\partial\hat{U}}{\partial\theta_{j}}\ket{\mathbf{x}_{p}}\bra{\mathbf{x}_{p}}\hat{U}^{\dagger}+\hat{U}\ket{\mathbf{x}_{p}}\bra{\mathbf{x}_{p}}\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{j}}\right)\,.$
(17)

Plugging Eqs.(16) and (17) in Eq.(7), we obtain the entries of the $M$. First
we compute $M_{\alpha_{p}\alpha_{q}}$, corresponding to the $\alpha_{p}$ and
$\alpha_{q}$ parameters:

$\begin{split}M_{\alpha_{p}\alpha_{q}}&=\mathrm{Tr}\left[{\frac{\partial\hat{\rho}^{\dagger}}{\partial\alpha_{p}}\frac{\partial\hat{\rho}}{\partial\alpha_{q}}}\right]\\\
&=\mathrm{Tr}\left[\hat{U}\ket{\mathbf{x}_{p}}\bra{\mathbf{x}_{p}}\hat{U}^{\dagger}\hat{U}\ket{\mathbf{x}_{q}}\bra{\mathbf{x}_{q}}\hat{U}^{\dagger}\right]\\\
&=\mathrm{Tr}\left[\hat{U}\ket{\mathbf{x}_{p}}\delta_{p,q}\bra{\mathbf{x}_{q}}\hat{U}^{\dagger}\right]\\\
&=\mathrm{Tr}\left[\ket{\mathbf{x}_{p}}\bra{\mathbf{x}_{q}}\hat{U}^{\dagger}\hat{U}\right]\delta_{p,q}\\\
&=\delta_{p,q}\delta_{p,q}=\delta_{p,q}\,,\end{split}$ (18)

where $\delta_{p,q}$ is the Kronecker delta. The remaining matrix elements
$M_{\alpha_{p}\theta_{j}}$ and $M_{\theta_{k}\theta_{j}}$ are then given by

$\begin{split}M_{\alpha_{p}\theta_{j}}&=\mathrm{Tr}\left[\frac{\partial\hat{\rho}^{\dagger}}{\partial\alpha_{p}}\frac{\partial\hat{\rho}}{\partial\theta_{j}}\right]\\\
&=\mathrm{Tr}\left[\hat{U}\ket{\mathbf{x}_{p}}\bra{\mathbf{x}_{p}}\hat{U}^{\dagger}\sum_{q}\alpha_{q}\left(\frac{\partial\hat{U}}{\partial\theta_{j}}\ket{\mathbf{x}_{q}}\bra{\mathbf{x}_{q}}\hat{U}^{\dagger}+\hat{U}\ket{\mathbf{x}_{q}}\bra{\mathbf{x}_{q}}\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{j}}\right)\right]\\\
&=\sum_{q}\alpha_{q}\delta_{p,q}\mathrm{Tr}\left[\frac{\partial\hat{U}}{\partial\theta_{j}}\ket{\mathbf{x}_{q}}\bra{\mathbf{x}_{p}}\hat{U}^{\dagger}+\hat{U}\ket{\mathbf{x}_{p}}\bra{\mathbf{x}_{q}}\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{j}}\right]\\\
&=\alpha_{p}\mathrm{Tr}\left[\ket{\mathbf{x}_{p}}\bra{\mathbf{x}_{p}}\hat{U}^{\dagger}\frac{\partial\hat{U}}{\partial\theta_{j}}+\ket{\mathbf{x}_{p}}\bra{\mathbf{x}_{p}}\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{j}}\hat{U}\right]\\\
&=\alpha_{p}\mathrm{Tr}\left[\ket{\mathbf{x}_{p}}\bra{\mathbf{x}_{p}}\left(\hat{U}^{\dagger}\frac{\partial\hat{U}}{\partial\theta_{j}}+\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{j}}\hat{U}\right)\right]\\\
&=\alpha_{p}\mathrm{Tr}\left[\ket{\mathbf{x}_{p}}\bra{\mathbf{x}_{p}}\frac{\partial\hat{U}^{\dagger}\hat{U}}{\partial\theta_{j}}\right]\\\
&=\alpha_{p}\mathrm{Tr}\left[\ket{\mathbf{x}_{p}}\bra{\mathbf{x}_{p}}\frac{\partial\mathbb{I}}{\partial\theta_{j}}\right]=0\,,\end{split}$
(19)

and

$\begin{split}M_{\theta_{k}\theta_{j}}&=\mathrm{Tr}\left[\frac{\partial\hat{\rho}^{\dagger}}{\partial\theta_{k}}\frac{\partial\hat{\rho}}{\partial\theta_{j}}\right]\\\
&=\mathrm{Tr}\left[\sum_{p}\alpha_{p}\left(\frac{\partial\hat{U}}{\partial\theta_{k}}\ket{\mathbf{x}_{p}}\bra{\mathbf{x}_{p}}\hat{U}^{\dagger}+\hat{U}\ket{\mathbf{x}_{p}}\bra{\mathbf{x}_{p}}\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{k}}\right)\sum_{q}\alpha_{q}\left(\frac{\partial\hat{U}}{\partial\theta_{j}}\ket{\mathbf{x}_{q}}\bra{\mathbf{x}_{q}}\hat{U}^{\dagger}+\hat{U}\ket{\mathbf{x}_{q}}\bra{\mathbf{x}_{q}}\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{j}}\right)\right]\\\
&=\mathrm{Tr}\left[\sum_{p,q}\alpha_{p}\alpha_{q}\left(\frac{\partial\hat{U}}{\partial\theta_{k}}\ket{\mathbf{x}_{p}}\bra{\mathbf{x}_{p}}\hat{U}^{\dagger}\frac{\partial\hat{U}}{\partial\theta_{j}}\ket{\mathbf{x}_{q}}{\bra{\mathbf{x}_{q}}}\hat{U}^{\dagger}+\frac{\partial\hat{U}}{\partial\theta_{k}}\ket{\mathbf{x}_{p}}\bra{\mathbf{x}_{q}}\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{j}}\delta_{p,q}+\delta_{p,q}\braket{\mathbf{x}_{p}}{\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{k}}\frac{\partial\hat{U}}{\partial\theta_{j}}}{\mathbf{x}_{q}}\right.\right.\\\
&+\left.\left.\hat{U}\ket{\mathbf{x}_{p}}\bra{\mathbf{x}_{p}}\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{k}}\hat{U}\ket{\mathbf{x}_{q}}\bra{\mathbf{x}_{q}}\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{j}}\right)\right]\\\
&=2\mathcal{R}e\left[\sum_{p}\alpha_{p}^{2}\bra{\mathbf{x}_{p}}\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{k}}\frac{\partial\hat{U}}{\partial\theta_{j}}\ket{\mathbf{x}_{p}}+\sum_{p,q}\alpha_{p}\alpha_{q}\bra{\mathbf{x}_{p}}\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{k}}\hat{U}\ket{\mathbf{x}_{q}}\bra{\mathbf{x}_{q}}\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{j}}\hat{U}\ket{\mathbf{x}_{p}}\right]\,.\end{split}$
(20)

The Liouvillian operator can be expanded, w.l.o.g. as
$\mathcal{L}=\sum_{r}c_{r}\hat{\sigma}_{r}^{1}\hat{\rho}\hat{\sigma}_{r}^{2}$,
with coefficients $c_{r}$ and Pauli strings $\hat{\sigma}_{r}^{1}$,
$\hat{\sigma}_{r}^{2}$. The components of $V$ can be then expressed as

$\begin{split}V_{\alpha_{p}}&=\mathrm{Tr}\left[\frac{\partial\hat{\rho}}{\partial\alpha_{p}}\mathcal{L}\left[\hat{\rho}\right]\right]\\\
&=\mathcal{R}e\left[\sum_{r,q}c_{r}\alpha_{q}\bra{\mathbf{x}_{p}}\hat{U}^{\dagger}\hat{\sigma}_{r}^{1}\hat{U}\ket{\mathbf{x}_{q}}\bra{\mathbf{x}_{q}}\hat{U}^{\dagger}\hat{\sigma}_{r}^{2}\hat{U}\ket{\mathbf{x}_{p}}\right]\,,\end{split}$
(21)

and

$\begin{split}V_{\theta_{k}}&=\mathrm{Tr}\left[\frac{\partial\hat{\rho}}{\partial\theta_{k}}\mathcal{L}\left[\hat{\rho}\right]\right]\\\
&=\mathcal{R}e\left[\sum_{rpq}c_{r}\alpha_{p}\alpha_{q}\left(\bra{\mathbf{x}_{p}}\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{k}}\hat{\sigma}_{r}^{1}\hat{U}\ket{\mathbf{x}_{q}}\bra{\mathbf{x}_{q}}\hat{U}^{\dagger}\hat{\sigma}_{r}^{2}\hat{U}\ket{\mathbf{x}_{p}}+\bra{\mathbf{x}_{p}}\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{k}}\hat{\sigma}_{r}^{2}\hat{U}\ket{\mathbf{x}_{q}}^{*}\bra{\mathbf{x}_{p}}\hat{U}^{\dagger}\hat{\sigma}_{r}^{1}\hat{U}\ket{\mathbf{x}_{q}}\right)\right]\,.\end{split}$
(22)

###  A.3 Equations of motion - Ansatz II

Proceeding similarly to section A.2, we derive the EOM of the parameters
$\bm{\beta}=(\bm{\alpha},\bm{\theta})$ by computing the elements $M_{kj}$ and
$V_{k}$ for Ansatz II
$\hat{\rho}(\bm{\beta)}=\sum_{p=1}^{R}\alpha_{p}\hat{U}(\bm{\theta}^{(p)})\ket{\mathbf{x}_{p}}\bra{\mathbf{x}_{p}}\hat{U}(\bm{\theta}^{(p)})^{\dagger}$,
with $\bm{\theta}=(\bm{\theta}^{(0)},\bm{\theta}^{(1)},...)$ be the set of
variational angular parameters. Again, we take for simplicity
$\hat{U}(\bm{\theta}^{(p)})=\hat{U}_{p}$. We begin by writing explicitly the
derivatives

$\frac{\partial\hat{\rho}}{\partial\alpha_{p}}=\hat{U}_{p}\ket{\mathbf{x}_{p}}\bra{\mathbf{x}_{p}}\hat{U}_{p}^{\dagger}\,,$
(23)
$\begin{split}\frac{\partial\hat{\rho}}{\partial\theta_{j}^{(q)}}=\alpha_{q}\left(\frac{\partial\hat{U}_{q}}{\partial\theta_{j}^{(q)}}\ket{\mathbf{x}_{q}}\bra{\mathbf{x}_{q}}\hat{U}_{q}^{\dagger}+\hat{U}_{q}\ket{\mathbf{x}_{q}}\bra{\mathbf{x}_{q}}\frac{\partial\hat{U}_{q}^{\dagger}}{\partial\theta_{j}^{(q)}}\right)\,.\end{split}$
(24)

Plugging these in Eq.(7), we obtain the elements $M_{\alpha_{p}\alpha_{q}}$
corresponding to the $\alpha_{p}$ and $\alpha_{q}$ parameters:

$\begin{split}M_{\alpha_{p}\alpha_{q}}&=|\bra{\mathbf{x}_{p}}\hat{U}_{p}^{\dagger}\hat{U}_{q}\ket{\mathbf{x}_{q}}|^{2}\,.\end{split}$
(25)

The remaining matrix elements $M_{\alpha_{p}\theta_{j}^{(q)}}$ and
$M_{\theta_{k}^{(p)}\theta_{j}^{(q)}}$ are simply

$\begin{split}M_{\alpha_{p}\theta_{j}^{(q)}}&=2\alpha_{p}\mathcal{R}e\left[\bra{\mathbf{x}_{p}}\hat{U}_{p}^{\dagger}\frac{\partial\hat{U}_{q}}{\partial\theta_{j}^{(q)}}\ket{\mathbf{x}_{q}}\bra{\mathbf{x}_{q}}\hat{U}_{q}^{\dagger}\hat{U}_{p}\ket{\mathbf{x}_{p}}\right]\,,\end{split}$
(26)

and

$\begin{split}M_{\theta_{k}^{(p)}\theta_{j}^{(q)}}&=2\alpha_{p}\alpha_{q}\mathcal{R}e\left[\bra{\mathbf{x}_{p}}\hat{U}_{p}^{\dagger}\hat{U}_{q}\ket{\mathbf{x}_{q}}\bra{\mathbf{x}_{q}}\frac{\partial\hat{U}_{q}^{\dagger}}{\partial\theta_{j}^{(q)}}\frac{\partial\hat{U}_{p}}{\partial\theta_{k}^{(p)}}\ket{\mathbf{x}_{p}}\right.\\\
&+\left.\bra{\mathbf{x}_{p}}\hat{U}_{p}^{\dagger}\frac{\partial\hat{U}_{q}}{\partial\theta_{j}^{(q)}}\ket{\mathbf{x}_{q}}\bra{\mathbf{x}_{q}}\hat{U}_{q}^{\dagger}\frac{\partial\hat{U}_{p}}{\partial\theta_{k}^{(p)}}\ket{\mathbf{x}_{p}}\right]\,.\end{split}$
(27)

Similarly, expanding the Liouvillian operator as
$\mathcal{L}[\hat{\rho}]=\sum_{r}c_{r}\hat{\sigma_{r}}^{1}\hat{\rho}\hat{\sigma_{r}}^{2}$,
the elements of the $V$ vector are given by

$\begin{split}V_{\alpha_{p}}=\mathcal{R}e\left[\sum_{r,q}c_{r}\alpha_{q}\bra{\mathbf{x}_{p}}\hat{U}_{p}^{\dagger}\hat{\sigma}_{r}^{1}\hat{U}_{q}\ket{\mathbf{x}_{q}}\bra{\mathbf{x}_{q}}\hat{U}_{q}^{\dagger}\hat{\sigma}_{r}^{2}\hat{U}_{p}\ket{\mathbf{x}_{p}}\right]\,,\end{split}$
(28)

and

$\begin{split}V_{\theta_{k}^{(p)}}&=\alpha_{p}\mathcal{R}e\left[\sum_{rq}c_{r}\alpha_{q}\left(\bra{\mathbf{x}_{p}}\frac{\partial\hat{U}_{p}^{\dagger}}{\partial\theta_{k}^{(p)}}\hat{\sigma}_{r}^{1}\hat{U}_{q}\ket{\mathbf{x}_{q}}\times\bra{\mathbf{x}_{q}}\hat{U}_{q}^{\dagger}\hat{\sigma}_{r}^{2}\hat{U}_{p}\ket{\mathbf{x}_{p}}\right.\right.\\\
&+\left.\left.\bra{\mathbf{x}_{p}}\frac{\partial\hat{U}_{p}^{\dagger}}{\partial\theta_{k}^{(p)}}\hat{\sigma}_{r}^{2}\hat{U}_{q}\ket{\mathbf{x}_{q}}^{*}\bra{\mathbf{x}_{p}}\hat{U}_{p}^{\dagger}\hat{\sigma}_{r}^{1}\hat{U}_{q}\ket{\mathbf{x}_{q}}\right)\right]\,.\end{split}$
(29)

##  Appendix B Truncated basis

When implementing LRQTE with a truncated basis, one must choose among the
$2^{n}$ possible basis states. In this work, we made this selection based on
the expansion of the Kraus operators $\hat{K}_{k}$ that govern the time
dynamics of the system. Let $\hat{\rho}$ be the density matrix of an open
quantum system. The time evolved state, given by the Lindblad master equation
Eq. (1), can be expressed as a Kraus map [45]

$\hat{\rho}(t+\delta
t)=\frac{1}{2K}\sum_{k=1}^{K}\left(\mathcal{\hat{U}}_{k}\hat{\rho}(t)\mathcal{\hat{U}}_{k}^{\dagger}+\mathcal{\hat{V}}_{k}\hat{\rho}(t)\mathcal{\hat{V}}_{k}^{\dagger}\right)+O(\delta
t^{2})\,,$ (30)

where $\\{\mathcal{\hat{U}}_{k}\\}_{k=1}^{K}$ and
$\\{\mathcal{\hat{V}}_{k}\\}_{k=1}^{K}$ are the Kraus operators of the
dynamics defined as

$\mathcal{\hat{U}}_{k}=e^{\delta t\hat{J}_{k}-i\sqrt{\gamma\delta
t}\hat{c}_{k}}\,,$ (31) $\mathcal{\hat{V}}_{k}=e^{\delta
t\hat{J}_{k}+i\sqrt{\gamma\delta t}\hat{c}_{k}}\,,$ (32)

with
$\hat{J}_{k}=-i\hat{H}-\frac{\gamma}{2}\left(\hat{c}_{k}^{2}-\hat{c}_{k}^{\dagger}\hat{c}_{k}\right)$.
Using the Zassenhaus approximation [65], we obtain

$\mathcal{\hat{U}}_{k}\hat{\rho}(t)\mathcal{\hat{U}}_{k}^{\dagger}+\mathcal{\hat{V}}_{k}\hat{\rho}(t)\mathcal{\hat{V}}_{k}^{\dagger}=e^{\delta
t\hat{J}_{k}}\hat{\Gamma}_{k}e^{\delta t\hat{J}_{k}^{\dagger}}+O(\delta
t^{2})\,,$ (33)

where $\hat{\Gamma}_{k}=\hat{\rho}(t)+\gamma\delta
t\left(\hat{c}_{k}\hat{\rho}(t)\hat{c}_{k}^{\dagger}-\hat{c}_{k}^{2}\hat{\rho}(t)\right)+h.c.$.
For the dissipative TFIM, with jump operators
$\hat{c}_{k}=\hat{\sigma}_{k}^{-}$, $\hat{c}_{k}^{2}=0$, we can further
simplify
$\hat{J}_{k}=-i\hat{H}+\frac{K}{2}\hat{\sigma}_{k}^{+}\hat{\sigma}_{k}^{-}$
and $\hat{\Gamma}_{k}=\hat{\rho}(t)+\gamma\delta
t\hat{\sigma}_{k}^{-}\hat{\rho}(t)\hat{\sigma}_{k}^{+}+h.c.$. Consider the
initial state $\hat{\rho}(0)=|1\rangle\langle 1|_{n}$. According to Eqs. (30)
and (33), the density matrix at a later time $\delta t$ is thus simply given
by

$\hat{\rho}(\delta t)=\frac{1}{K}\sum_{k=1}^{K}\left(e^{-i\delta
t\hat{H}}\hat{\rho}(0)e^{i\delta t\hat{H}^{\dagger}}\right)+O(\delta
t^{2})\,,$ (34)

where we have used $\hat{\sigma}_{k}^{-}|1\rangle\langle
1|_{n}\hat{\sigma}_{k}^{+}=0$. Eq. (34) describes the Hamiltonian dynamics of
a quantum system isolated from its environment. In this picture, the time-
evolved state of the system is approximately

$\begin{split}|\psi(t+\delta t)\rangle&\approx\left(\mathbb{I}-i\delta
t(\hat{H}_{x}+\hat{H}_{zz})\right)|\psi(t)\rangle\\\
&\in\textrm{span}\left(\\{|\psi(t)\rangle,\hat{\sigma}_{1}^{x}|\psi(t)\rangle,\hat{\sigma}_{2}^{x}|\psi(t)\rangle,...,\hat{\sigma}_{n}^{x}|\psi(t)\rangle\\}\right)\,,\end{split}$
(35)

where we have used that the diagonal component of the Hamiltonian
$\hat{H}_{zz}$ is diagonal in the computational basis, acting trivially on
$|\psi\rangle$. Therefore, a natural choice of basis for $|\psi(t+\delta
t)\rangle$ is given by the quantum states at Hamming distance 1 from
$|\psi(t)\rangle$. Analogously, $|\psi(t+2\delta t)\rangle$ will be in a
subspace generated by the states at Hamming distance $\leq 2$ which can be
used for an extension of the basis. In the following, we investigate how the
choice of basis affects the performance of Ansätze I and II. Fig.8 shows the
infidelity of a 4-spin (2D) TFIM system. The simulation parameters are
$J_{z}=\gamma=1$, $h=0.5$, $dt=5\times 10^{-3}$, $N_{l}=4$ and the rank
$R=\\{5,16\\}$. Furthermore, the states $|\mathbf{x}_{p}\rangle$ of the
truncated rank simulations are taken from a basis $\mathcal{A}$

$\begin{split}\mathcal{A}=\\{|1111\rangle,|1110\rangle,|1100\rangle,|1101\rangle,|1011\rangle\\}\,,\end{split}$

and a basis $\mathcal{B}$, built from the states at Hamming distance 1 from
the initial state,

$\begin{split}\mathcal{B}=\\{|1111\rangle,|1110\rangle,|1101\rangle,|1011\rangle,|0111\rangle\\}\,.\end{split}$

As expected, the basis $\mathcal{B}$ leads to an overall lower infidelity. In
particular, the low-rank results computed with this basis for Ansatz I are
significantly better than the ones obtained with basis $\mathcal{A}$. On the
other hand, the choice of the initial basis doesn’t seem to affect the
performance of the Ansatz II. Naturally, this is due to the fact that, in this
representation, each basis vector is individually rotated according to a
unitary $\hat{U}_{i}$. This feature significantly increases the flexibility,
and therefore also the expressivity of the variational Ansatz. In contrast,
LRQTE with Ansatz I optimizes a single parameterized unitary operator that is
applied to all basis vectors, and, therefore, the initial choice of the basis
states is significantly more important.

Figure 8: Influence of the initial basis choice on the performance of the
LRQTE algorithm for Ansatz I and II.

##  Appendix C Resource requirements

In this section, we analyse the resource requirements of the proposed LRQTE,
both for Ansatz I and II. In this analysis, we focus on the quantum cost of
each time step of the algorithm which we express in terms of the number of
expectation values to compute. Moreover, by deriving the number of expectation
values that we must compute, we are better informed about the error in the
estimation of the variational parameters when estimated in a real quantum
device, since the error in the computation of the expectation values
propagates into the elements of $M$ and $V$ in Eq. (7). In tables 1 and 2 we
present the number of expectation values to be computed in one iteration of
the algorithm, which we express in terms of the rank $R$ and the number of
layers $N_{l}$ for each variational quantum circuit. Note that for Ansatz I
the total number of angular parameters $N_{I,\theta}\propto N_{l}$ and for
Ansatz II, $N_{II,\theta}\propto RN_{l}$. We conclude that the computational
cost of implementing Ansatz I and II is, respectively,
$O(R^{2}N_{l}+RN_{l}^{2})$ and $O(R^{2}N_{l}^{2})$. As a final remark, notice
that while the full-rank simulation with rank $R=2^{n}$ would lead to an
exponential number of terms to evaluate, the proposed low-rank Ansätze require
only $R=O(\mbox{poly}(n))$, keeping the computational cost polynomial with the
size of the system.

| Number of Expectation Values
---|---
$M_{\alpha_{p},\alpha_{q}}$ | $0$
$M_{\alpha_{p},\theta_{j}}$ | $0$
$M_{\theta_{k},\theta_{j}}$ | $O(R^{2})$
$\bm{M}$ | $\bm{O(R^{2}N_{l}+RN_{l}^{2})}$
$V_{\alpha_{p}}$ | $O(R)$
$V_{\theta_{k}}$ | $O(R^{2})$
$\bm{V}$ | $\bm{O(R^{2}N_{l})}$
Table 1: Number of expectation values needed to estimate $M$ and $V$ for Ansatz I. | Number of Expectation Values
---|---
$M_{\alpha_{p},\alpha_{q}}$ | $O(1)$
$M_{\alpha_{p},\theta_{j}^{(q)}}$ | $O(1)$
$M_{\theta_{k}^{(p)},\theta_{j}^{(q)}}$ | $O(1)$
$\bm{M}$ | $\bm{O(R^{2}N_{l}^{2})}$
$V_{\alpha_{p}}$ | $O(R)$
$V_{\theta_{k}^{(p)}}$ | $O(R)$
$\bm{V}$ | $\bm{O(R^{2}N_{l})}$
Table 2: Number of expectation values needed to estimate $M$ and $V$ for
Ansatz II.

##  Appendix D Quantum circuits

The numerical values of $M$ and $V$ in Eq. (7) are obtained by evaluating
expectation values of quantum circuits. In this section, we describe how we
compute these on a quantum device for Ansatz I. The same discussion applies to
Ansatz II.

The variational unitary operator can be expressed as a sequence of
parameterized circuits

$\hat{U}(\bm{\theta})=\hat{U}_{N}(\theta_{N})\hat{U}_{N-1}(\theta_{N-1})...\hat{U}_{1}(\theta_{1})\,,$
(36)

where
$\hat{U}_{k}({\theta_{k}})=\hat{W}_{k}e^{-i\frac{\theta_{k}}{2}\hat{H}_{k}}$,
$\hat{W}_{k}$s are constant operators and $\hat{H}_{k}$ are Hermitian
operators which can be expressed as the sum of Pauli strings

$\hat{H}_{k}=\sum_{s}g_{k,s}\hat{\sigma}_{k,s}\,,$ (37)

with $g_{k,s}\in\mathbb{R}$. It follows that

$\displaystyle\frac{\partial\hat{U}}{\partial\theta_{k}}|\mathbf{x}_{p}\rangle$
$\displaystyle=\hat{U}_{N}(\theta_{N})...\frac{\partial\hat{U}_{k}(\theta_{k})}{\partial\theta_{k}}...\hat{U}_{1}(\theta_{1})|\mathbf{x}_{p}\rangle$
(38)
$\displaystyle=\sum_{s}-\frac{i}{2}g_{k,s}\widetilde{U}_{k,s}\ket{\mathbf{x}_{p}}\,,$

with

$\displaystyle\widetilde{U}_{k,s}$
$\displaystyle=\hat{U}_{N:k}\hat{\sigma}_{k,s}\hat{U}_{k-1:1}\,,$ (39)
$\displaystyle\hat{U}_{N:k}$
$\displaystyle=\hat{U}_{N}(\theta_{N})...\hat{U}_{k}(\theta_{k})\,,$
$\displaystyle\hat{U}_{k-1:1}$
$\displaystyle=\hat{U}_{k-1}(\theta_{k-1})...\hat{U}_{1}(\theta_{1})\,.$

From Eqs. (20), (21) and (22) we conclude that we need to compute terms of the
form$\\{\braket{\mathbf{x}_{p}}{\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{k}}\frac{\partial\hat{U}}{\partial\theta_{j}}}{\mathbf{x}_{p}}$,
$\braket{\mathbf{x}_{p}}{\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{k}}\hat{U}}{\mathbf{x}_{q}}$,
$\braket{\mathbf{x}_{p}}{\hat{U}^{\dagger}\hat{\sigma}_{r}\hat{U}}{\mathbf{x}_{q}}$,
$\braket{\mathbf{x}_{p}}{\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{k}}\hat{\sigma}_{r}\hat{U}}{\mathbf{x}_{q}}\\}$.
Take
$\bra{\mathbf{x}_{p}}\frac{\partial\hat{U}^{{\dagger}}}{\partial\theta_{k}}\hat{U}\ket{\mathbf{x}_{q}}$
as an example. It can be expanded as

$\displaystyle\bra{\mathbf{x}_{p}}\frac{\partial\hat{U}^{{\dagger}}}{\partial\theta_{k}}\hat{U}\ket{\mathbf{x}_{q}}$
$\displaystyle=\frac{i}{2}\sum_{s}g_{k,s}\bra{\mathbf{x}_{p}}\widetilde{U}_{k,s}^{{\dagger}}\hat{U}\ket{\mathbf{x}_{q}}$
(40)
$\displaystyle=\frac{i}{2}\sum_{s}g_{k,s}\bra{\mathbf{x}_{p}}\hat{U}_{k-1:1}^{{\dagger}}\hat{\sigma}_{k,s}\hat{U}_{k-1:1}\ket{\mathbf{x}_{q}}$
$\displaystyle=\frac{i}{2}\sum_{s}g_{k,s}\bra{\mathbf{x}_{q}}\hat{X}_{p-q}\hat{U}_{k-1:1}^{{\dagger}}\hat{\sigma}_{k,s}\hat{U}_{k-1:1}\ket{\mathbf{x}_{q}}$
$\displaystyle=\frac{1}{2}\sum_{s}g_{k,s}(\mathcal{I}m[\bra{\mathbf{x}_{q}}\hat{U}_{k-1:1}^{{\dagger}}\hat{\sigma}_{k,s}\hat{U}_{k-1:1}\hat{X}_{p-q}\ket{\mathbf{x}_{q}}]$
$\displaystyle+i\mathcal{R}e[\bra{\mathbf{x}_{q}}\hat{U}_{k-1:1}^{{\dagger}}\hat{\sigma}_{k,s}\hat{U}_{k-1:1}\hat{X}_{p-q}\ket{\mathbf{x}_{q}}])\,,$

where $\hat{X}_{p-q}$ is the permutation gate that satisfies
$\ket{\mathbf{x}_{p}}=\hat{X}_{p-q}\ket{\mathbf{x}_{q}}$. Therefore, we need
to evaluate terms of the form
$\langle\psi|\hat{V}^{{\dagger}}\hat{\sigma}_{2}\hat{V}\hat{\sigma}_{1}|\psi\rangle$,
where $\hat{V}$ is a known unitary operator, and $\hat{\sigma}_{2}$,
$\hat{\sigma}_{1}$ are Pauli strings. In Eq. (40),
$\ket{\psi}=\ket{\mathbf{x}_{q}}$, $\hat{V}=\hat{U}_{k-1:1}$,
$\hat{\sigma}_{1}=\hat{X}_{p-q}$, $\hat{\sigma}_{2}=\hat{\sigma}_{k,s}$. We
consider two strategies to evaluate the value above. The Hadamard test [48]
requires one additional (ancilla) qubit (see Fig.9) and it outputs the real
and imaginary parts of
$\langle\psi|\hat{V}^{{\dagger}}\hat{\sigma}_{2}\hat{V}\hat{\sigma}_{1}|\psi\rangle$,
by setting the angle of rotation gate around the $z$-axis to 0 or $-\pi/2$,
respectively.

Figure 9: Hadamard test circuit for estimating
$\langle\psi|\hat{V}^{{\dagger}}\hat{\sigma}_{2}\hat{V}\hat{\sigma}_{1}|\psi\rangle$.
The first register with input $\ket{0}$ is an ancilla qubit. Setting the angle
$\theta_{z}$ of the $R_{z}=e^{-i\frac{\theta_{z}}{2}\hat{\sigma}_{z}}$ gate to
0 or $-\pi/2$ outputs the real or imaginary part of
$\langle\psi|\hat{V}^{{\dagger}}\hat{\sigma}_{2}\hat{V}\hat{\sigma}_{1}|\psi\rangle$,
respectively. The real (and imaginary) part of the expectation value is
calculated as $2P(0)-1$, where $P(0)$ is the probability that we measure the
state $|0\rangle$. Figure 10: Ancilla-free quantum circuit that evaluates
$\mathcal{R}e\left[\langle\psi|\hat{V}^{{\dagger}}_{2}\hat{V}_{1}|\psi\rangle\right]$.
The first measurement prepares the state
$\hat{P}_{1,\pm}|\psi\rangle/\lVert\hat{P}_{1,\pm}|\psi\rangle\rVert$ and the
second measurement estimates $\langle\hat{\sigma}_{2}\rangle$. Figure 11:
Ancilla-free quantum circuit that evaluates
$\mathcal{I}m\left[\langle\psi|\hat{V}^{{\dagger}}_{2}\hat{V}_{1}|\psi\rangle\right]$.
Two circuits are required, one with the multi-qubit gate
$e^{i\pi\frac{\hat{\sigma}_{1}}{4}}$ and the other with
$e^{-i\pi\frac{\hat{\sigma}_{1}}{4}}$.

When implementing the Hadamard test on a superconducting quantum computer, the
limited geometry of current quantum devices forces the addition of faulty two
qubit SWAP gates in order to implement operations between qubits that are not
physically connected. For this reason, the controlled gates required by this
method increase considerably the probability of error.

An alternative to the Hadamard test that requires no additional qubits is
proposed in [66]. To illustrate this method, we rewrite the term
$\langle\psi|\hat{V}^{{\dagger}}\hat{\sigma}_{2}\hat{V}\hat{\sigma}_{1}|\psi\rangle$
as

$\mathcal{R}e\left[\langle\psi|\hat{V}^{{\dagger}}\hat{\sigma}_{2}\hat{V}\hat{\sigma}_{1}|\psi\rangle\right]=\langle\psi|\hat{P}_{1,+}\hat{V}^{{\dagger}}\hat{\sigma}_{2}\hat{V}\hat{P}_{1,+}|\psi\rangle-\langle\psi|\hat{P}_{1,-}\hat{V}^{{\dagger}}\hat{\sigma}_{2}\hat{V}\hat{P}_{1,-}|\psi\rangle\,,$
(41)

where $\hat{P}_{1,\pm}=(1\pm\hat{\sigma}_{1})/2$ is the projection operator
that measures the Pauli string $\hat{\sigma}_{1}$. Instead of implementing
this multi-qubit measurements of $\hat{\sigma}_{1}$ and $\hat{\sigma}_{2}$, we
replace it by a single-qubit measurement by applying the appropriate basis
change operators $\hat{V}_{1}$ and $\hat{V}_{2}$. The quantum circuit for
evaluating
$\mathcal{R}e[\langle\psi|\hat{V}^{{\dagger}}\hat{\sigma}_{2}\hat{V}\hat{\sigma}_{1}|\psi\rangle]$
is provided in Fig.10. Notice that after the projective measurement, the
resulting normalized state is
$\hat{P}_{1,\pm}|\psi\rangle/\lVert\hat{P}_{1,\pm}|\psi\rangle\rVert$, where
$\lVert\hat{P}_{1,\pm}|\psi\rangle\rVert$ is estimated from the probability
distribution of the projective measurement. As for the imaginary part, the
ancilla-free quantum circuit is derived from the following formula

$\displaystyle\mathcal{I}m[\langle\psi|\hat{V}^{{\dagger}}\hat{\sigma}_{2}\hat{V}\hat{\sigma}_{1}|\psi\rangle]=\frac{1}{2}(\langle\psi|e^{-i\pi\frac{\hat{\sigma}_{1}}{4}}\hat{V}^{{\dagger}}\hat{\sigma}_{2}\hat{V}e^{i\pi\frac{\hat{\sigma}_{1}}{4}}|\psi\rangle$
(42)
$\displaystyle-\langle\psi|e^{i\pi\frac{\hat{\sigma}_{1}}{4}}\hat{V}^{{\dagger}}\hat{\sigma}_{2}\hat{V}e^{-i\pi\frac{\hat{\sigma}_{1}}{4}}|\psi\rangle)\,.$

Therefore, the imaginary part can be estimated with two quantum circuits (see
Fig.11), which measure the expectation value of $\hat{\sigma}_{2}$ over the
states $|\Psi_{\pm}\rangle=\hat{V}e^{\pm
i\pi\frac{\hat{\sigma}_{1}}{4}}|\psi\rangle$.

Figure 12: Hadamard test circuit for computing
$\braket{\psi}{\hat{W}^{\dagger}\hat{V}^{\dagger}\hat{\sigma}_{3}\hat{V}\hat{\sigma}_{2}\hat{W}\hat{\sigma}_{1}}{\psi}$.
Setting the angle $\theta_{z}$ of the $\hat{R}_{z}$ gate as 0 or $-\pi/2$
measures the real and imaginary part, respectively. The value is estimated as
$2P(0)-1$, where $P(0)$ is the probability that the ancilla qubit is measured
to be 0.

Similarly, we can show that the circuits displayed in Figs.9, 10, and 11 can
be used to compute
$\braket{\mathbf{x}_{p}}{\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{k}}\frac{\partial\hat{U}}{\partial\theta_{j}}}{\mathbf{x}_{p}}$,
$\braket{\mathbf{x}_{p}}{\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{k}}\hat{U}}{\mathbf{x}_{q}}$,
and
$\braket{\mathbf{x}_{p}}{\hat{U}^{\dagger}\hat{\sigma}_{r}\hat{U}}{\mathbf{x}_{q}}$.
In order to evaluate
$\braket{\mathbf{x}_{p}}{\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{k}}\hat{\sigma}_{r}\hat{U}}{\mathbf{x}_{q}}$
with an optimal circuit depth, 3 control gates are required in the Hadamard
test, as shown in Fig.12. The expectation value
$\bra{\mathbf{x}_{p}}\frac{\partial\hat{U}^{{\dagger}}}{\partial\theta_{k}}\hat{\sigma}_{r}\hat{U}\ket{\mathbf{x}_{q}}$
can be written as

$\displaystyle\bra{\mathbf{x}_{p}}\frac{\partial\hat{U}^{{\dagger}}}{\partial\theta_{k}}\hat{\sigma}_{r}\hat{U}\ket{\mathbf{x}_{q}}=\frac{i}{2}\sum_{s}g_{k,s}\bra{\mathbf{x}_{p}}\widetilde{U}_{k,s}^{{\dagger}}\hat{U}\ket{\mathbf{x}_{q}}$
(43)
$\displaystyle=\frac{i}{2}\sum_{s}g_{k,s}\bra{\mathbf{x}_{p}}\hat{U}_{k-1:1}^{{\dagger}}\hat{\sigma}_{k,s}\hat{U}_{N:k}^{{\dagger}}\hat{\sigma}_{r}\hat{U}_{N:k}\hat{U}_{k-1:1}\ket{\mathbf{x}_{q}}$
$\displaystyle=\frac{i}{2}\sum_{s}g_{k,s}\bra{\mathbf{x}_{q}}\hat{X}_{p-q}\hat{U}_{k-1:1}^{{\dagger}}\hat{\sigma}_{k,s}\hat{U}_{N:k}^{{\dagger}}\hat{\sigma}_{r}\hat{U}_{N:k}\hat{U}_{k-1:1}\ket{\mathbf{x}_{q}}$
$\displaystyle=\frac{i}{2}\sum_{s}g_{k,s}\bra{\mathbf{x}_{q}}\hat{U}_{k-1:1}^{{\dagger}}\hat{U}_{N:k}^{{\dagger}}\hat{\sigma}_{r}\hat{U}_{N:k}\hat{\sigma}_{k,s}\hat{U}_{k-1:1}\hat{X}_{p-q}\ket{\mathbf{x}_{q}}^{*}\,,$

which is of the form
$\braket{\psi}{\hat{W}^{\dagger}\hat{V}^{\dagger}\hat{\sigma}_{3}\hat{V}\hat{\sigma}_{2}\hat{W}\hat{\sigma}_{1}}{\psi}$,
where $\ket{\psi}=\ket{\mathbf{x}_{q}}$, $\hat{\sigma}_{1}=\hat{X}_{p-q}$,
$\hat{\sigma}_{2}=\hat{\sigma}_{k,s}$, $\hat{\sigma}_{3}=\hat{\sigma}_{r}$,
$\hat{V}=\hat{U}_{N:k}$, and $\hat{W}=\hat{U}_{k-1:1}$. Evaluating
$\braket{\psi}{\hat{W}^{\dagger}\hat{V}^{\dagger}\hat{\sigma}_{3}\hat{V}\hat{\sigma}_{2}\hat{W}\hat{\sigma}_{1}}{\psi}$
also requires more circuits than
$\langle\psi|\hat{V}^{{\dagger}}\hat{\sigma}_{2}\hat{V}\hat{\sigma}_{1}|\psi\rangle$
(more details in [66]). In summary, while the Hadamard test circuit requires
one additional qubit, and potentially many additional noisy SWAP gates, the
ancilla-free method requires more circuits. The choice between the two schemes
then depends on several factors, namely the exact architecture of the
superconducting device, the error rates, and the system size. In this work, we
opted for an hybrid approach of these two methods, in which we compute the
terms
$\braket{\mathbf{x}_{p}}{\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{k}}\frac{\partial\hat{U}}{\partial\theta_{j}}}{\mathbf{x}_{p}}$,
$\braket{\mathbf{x}_{p}}{\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{k}}\hat{U}}{\mathbf{x}_{q}}$,
$\braket{\mathbf{x}_{p}}{\hat{U}^{\dagger}\hat{\sigma}_{r}\hat{U}}{\mathbf{x}_{q}}$
with ancilla-free circuits and
$\braket{\mathbf{x}_{p}}{\frac{\partial\hat{U}^{\dagger}}{\partial\theta_{k}}\hat{\sigma}_{r}\hat{U}}{\mathbf{x}_{q}}$
with the Hadamard test. In comparison to the method that computes all
expectation values with the Hadamard test, this scheme improves the accuracy
of the results, while keeping the computational resources manageable.

##  Appendix E Error analysis

###  E.1 Density leakage in the low-rank approximation

In the LRQTE, the trace of the hybrid-quantum classical Ansatz which
represents the variational density matrix, is simply $\sum_{p}{\alpha_{p}}$.
In this section, we show that, while the trace is preserved for the full rank
simulation, it decreases in the low-rank dynamics, due to the incomplete
representation space that it can generate. This analysis is presented here for
Ansatz I, but can be generalized to any low-rank scheme [49]. From Eqs. (18),
(19) and (21), we obtain

$\dot{\alpha}_{p}=\mathrm{Tr}\left[\hat{U}\ket{\mathbf{x}_{p}}\bra{\mathbf{x}_{p}}\hat{U}^{\dagger}\mathcal{L}\left[\hat{\rho}\right]\right]=\braket{\mathbf{x}_{p}}{\hat{U}^{\dagger}\mathcal{L}[\hat{\rho}]\hat{U}}{\mathbf{x}_{p}}\,.$
(44)

Therefore, the time derivative of the trace is given by

$\mathrm{Tr}\left[\dot{\hat{\rho}}(\bm{\beta})\right]=\sum_{p}\dot{\alpha_{p}}=\sum_{p}\braket{\mathbf{x}_{p}}{\hat{U}^{\dagger}\mathcal{L}[\hat{\rho}]\hat{U}}{\mathbf{x}_{p}}\,,$
(45)

where $|\mathbf{x}_{p}\rangle$ is a computational basis state. Using that for
the full rank Ansatz $\\{|\mathbf{x}_{p}\rangle\\}_{p=1}^{R=2^{n}}$ is a
complete basis set that satisfies
$\sum_{p}|\mathbf{x}_{p}\rangle\langle\mathbf{x}_{p}|=\mathbb{I}$, it follows
that

$\displaystyle\mathrm{Tr}\left[\dot{\hat{\rho}}_{full}\right]$
$\displaystyle=\sum_{p}\braket{\mathbf{x}_{p}}{\hat{U}^{\dagger}\mathcal{L}[\hat{\rho}_{full}]\hat{U}}{\mathbf{x}_{p}}$
(46)
$\displaystyle=\mathrm{Tr}\left[\hat{U}^{{\dagger}}\mathcal{L}[\hat{\rho}_{full}]\hat{U}\right]$
$\displaystyle=\mathrm{Tr}[\mathcal{L}[\hat{\rho}_{full}]]=0\,,$

where we’ve used the trace-preserving propriety of the Liouvillian super
operator. Therefore, we conclude that the trace is preserved for the full rank
variational time evolution.

Let the computational basis states that are included and excluded in the low-
rank Ansatz with subscripts $q$ and $q^{\prime}$, respectively. We can
decompose Eq.(46) into two sums over the $|\mathbf{x}_{q}\rangle$ and the
$|\mathbf{x}_{q^{\prime}}\rangle$ states:

$\sum_{q}\braket{\mathbf{x}_{q}}{\hat{U}^{\dagger}\mathcal{L}[\hat{\rho}]\hat{U}}{\mathbf{x}_{q}}+\sum_{q^{\prime}}\braket{\mathbf{x}_{q^{\prime}}}{\hat{U}^{\dagger}\mathcal{L}[\hat{\rho}]\hat{U}}{\mathbf{x}_{q^{\prime}}}=0\,.$
(47)

Noticing that the sum over the $|\mathbf{x}_{q}\rangle$ states is the low-rank
density matrix, we can conclude that

$\mathrm{Tr}\left[\dot{\hat{\rho}}_{low}\right]=\sum_{q}\braket{\mathbf{x}_{q}}{\hat{U}^{\dagger}\mathcal{L}[\hat{\rho}_{low}]\hat{U}}{\mathbf{x}_{q}}=-\sum_{q^{\prime}}\braket{\mathbf{x}_{q^{\prime}}}{\hat{U}^{\dagger}\mathcal{L}[\hat{\rho}_{low}]\hat{U}}{\mathbf{x}_{q^{\prime}}}\,.$
(48)

Since
$\hat{\rho}_{low}\hat{U}|\mathbf{x}_{q^{\prime}}\rangle=\sum_{q}\alpha_{q}\hat{U}\ket{\mathbf{x}_{q}}\bra{\mathbf{x}_{q}}\hat{U}^{\dagger}\hat{U}|\mathbf{x}_{q^{\prime}}\rangle=0$,
we obtain

$\mathrm{Tr}\left[\dot{\hat{\rho}}_{low}\right]=-\sum_{q^{\prime}}\braket{\mathbf{x}_{q^{\prime}}}{\hat{U}^{\dagger}\mathcal{D}[\hat{\rho}_{low}]\hat{U}}{\mathbf{x}_{q^{\prime}}}\,,$
(49)

with
$\mathcal{D}[\hat{\rho}]=\sum_{k}\gamma_{k}\hat{c}_{k}\hat{\rho}\hat{c}_{k}^{\dagger}$.
It can be easily proved that$\mathrm{Tr}\left[\dot{\rho}_{low}\right]\leq 0$.
We can then conclude that, whereas the algorithm preserves the trace of the
full-rank variational density matrix, the same cannot be guaranteed for the
low-rank simulation.

###  E.2 Error bounds

The calculation of the infidelity requires the knowledge of the exact time
evolved state, whose computation is not tractable for large-scale simulations.
Therefore, it is useful to derive a post-processing error bound based on the
data from the LRQTE, for which the exact time evolution is not required. We
follow the method introduced in [67] to construct this error bound. Let the
density matrix of LRQTE be $\hat{\rho}(\bm{\beta}(t))$ and the exact time-
evolved state $\hat{\rho}^{*}(t)$. The $L^{2}$-distance between
$\hat{\rho}(\bm{\beta}(t))$ and $\hat{\rho}^{*}(t)$ is

$\lVert\hat{\rho}(\bm{\beta}(t))-\hat{\rho}^{*}(t)\rVert_{2}=\sqrt{\mathrm{Tr}[\hat{\rho}(\bm{\beta})^{2}]+\mathrm{Tr}[\hat{\rho}^{*2}]-2\mathrm{Tr}[\hat{\rho}(\bm{\beta})\hat{\rho}^{*}]}\,.$
(50)

Since the Mclachlan’s variational principle is obtained from the minimization
of the $L^{2}$-distance
$\lVert\sum_{j}\frac{\partial\hat{\rho}(\bm{\beta})}{\partial\beta_{j}}\dot{\beta_{j}}-\mathcal{L}\left[\hat{\rho}(\bm{\beta})\right]\rVert$,
we can define a cost function for the LRQTE at each time step given by Eq.
(14) and expand it as

$\begin{split}&C_{Mcl}(\dot{\bm{\beta}})=\sum_{j,k}\mathrm{Tr}\left[\frac{\partial\hat{\rho}(\bm{\beta})}{\partial\beta_{j}}\frac{\partial\hat{\rho}(\bm{\beta})}{\partial\beta_{k}}\right]\dot{\beta}_{j}\dot{\beta}_{k}\\\
&-2\sum_{j}Tr\left[\frac{\partial\hat{\rho}(\bm{\beta})}{\partial\beta_{j}}\mathcal{L}[\hat{\rho}(\bm{\beta})]\right]\dot{\beta}_{j}+\mathrm{Tr}\left[\mathcal{L}[\hat{\rho}(\bm{\beta})]^{2}\right]\\\
&=\sum_{j,k}M_{j,k}\dot{\beta}_{j}\dot{\beta}_{k}-2\sum_{j}V_{j}\dot{\beta}_{j}+\mathrm{Tr}\left[\mathcal{L}[\hat{\rho}(\bm{\beta})]^{2}\right]\,.\end{split}$
(51)

Note that, by minimizing $C_{Mcl}$ with respect to $\dot{\bm{\beta}}$, we
recover Eq.(7). Moreover, we can rewrite the cost function at time $t$ as a
function of the variables $(\bm{\beta}(t),\bm{\dot{\beta}}(t))$ by considering
$M_{j,k}$ and $V_{k}$ as functions of $\bm{\beta}(t)$

$\begin{split}&C(\bm{\beta}(t),\bm{\dot{\beta}}(t))=\sum_{j,k}M_{j,k}(\bm{\beta}(t))\dot{\beta}_{j}(t)\dot{\beta}_{k}(t)\\\
&-2\sum_{j}V_{j}(\bm{\beta}(t))\dot{\beta}_{j}(t)+\mathrm{Tr}\left[\mathcal{L}[\hat{\rho}(\bm{\beta}(t))]^{2}\right]\,,\end{split}$
(52)

where the values of $M_{j,k}(\bm{\beta}(t))$, $V_{j}(\bm{\beta}(t))$ and
$\dot{\bm{\beta}}(t)$ are directly obtained from the LRQTE algorithm, while
$\mathrm{Tr}\left[\mathcal{L}[\hat{\rho}(\bm{\beta})(t))]^{2}\right]$ needs
additional measurements. We define the posterior error bound of LRQTE as

$E_{p}(T)=\int_{0}^{T}\sqrt{C(\bm{\beta}(t),\bm{\dot{\beta}}(t))}dt\,.$ (53)

It’s easy to show that
$\lVert\hat{\rho}(\bm{\beta}(T)))-\hat{\rho}^{*}(T)\rVert_{2}\leq E_{p}(T)$,
and therefore, $E_{p}(T)$ provides an upper bound of the $L^{2}$-distance
between the LRQTE and exact time evolution states. When computing the value of
$E_{p}$ from the simulated variational parameters, the integral in Eq.(53) can
be approximated by a sum over the discretized cost function, given by

$\begin{split}C_{i}(\bm{\beta}(t_{i}),\bm{\beta}(t_{i+1}))&=\sum_{j,k}M_{j,k}(\bm{\beta}(t_{i}))\frac{\beta_{j}(t_{i+1})-\beta_{j}(t_{i})}{t_{i+1}-t_{i}}\frac{\beta_{k}(t_{i+1})-\beta_{k}(t_{i})}{t_{i+1}-t_{i}}\\\
&-2\sum_{j}V_{j}(\bm{\beta}(t_{i}))\frac{\beta_{j}(t_{i+1})-\beta_{j}(t_{i})}{t_{i+1}-t_{i}}+\mathrm{Tr}\left[\mathcal{L}[\hat{\rho}(\bm{\beta}(t_{i}))]^{2}\right]\,.\end{split}$
(54)

Then the posterior error bound $E_{p}$ is

$\begin{split}E_{p}({t_{0},t_{1},...t_{N}},{\bm{\beta}(t_{0}),\bm{\beta}(t_{1}),...\bm{\beta}(t_{N})})=\sum_{i=0}^{N-1}\sqrt{C_{i}\times(t_{i+1}-t_{i})}\,.\end{split}$
(55)

Fig.13 shows the $E_{p}$ of the LRQTE computed on different simulators and the
$L^{2}$-distance between the state given by LRQTE algorithm and the exact one.
Notice that, although the estimation of the variational parameters is
performed on different simulators, the calculation of the posterior error
bound and the $L^{2}$-distance are all obtained from a statevector simulation,
which gives an accurate estimation of the cost function. Fig.13, obtained by
simulating a 2-site TFIM with $J_{z}=\gamma=1$, $h=0.5$, $dt=1\times 10^{-2}$
and $R=4$, verifies that the posterior error defined in Eq.(55) upper bounds
the $L^{2}$-distance, and can be used as an estimator of the accuracy of the
LRQTE.

Figure 13: Posterior error bound for the 2-site TFIM LRQTE simulation and the
$L^{2}$-distance between LRQTE and the exact time evolution.

##  Appendix F Numerical regularization schemes

When implementing the LRQTE, the linear system Eq.(7) could be ill-conditioned
during the simulation, i.e., when $M$ is singular or it’s eigenvalues
extremely small, inverting it may cause numerical instability when solving the
linear system, leading to errors in the estimation of the variational
parameters. In this work, we address this issue by adopting different kinds of
regularization schemes for both the statevector and the Qiskit simulations.
For the statevector simulations, two regularization methods were implemented.
One consisted in truncating the smallest eigenvalues and solving the reduced
linear system (which is more stable), and the second, in rescaling the
eigenvalues of the pseudoinverse of $M$ [68]. The latter can be summarized as
follows: (1) first, we diagonalize the matrix $M=U\Sigma U^{\dagger}$, (2)
then, find the eigenvalues $\sigma_{\mu}^{2}$ of
$\Sigma=\textrm{diag}(\sigma_{1}^{2},...,\sigma_{P}^{2})$, (3) finally, define
the pseudoinverse $\tilde{M}^{-1}=U\tilde{\Sigma}^{-1}U^{\dagger}$, with

$\tilde{\Sigma}_{\mu,\nu}^{-1}=\frac{1/\hat{\sigma}_{\mu}^{2}}{1+(\lambda^{2}/\hat{\sigma}_{\mu}^{2})^{6}}\delta_{\mu,\nu}\,,$
(56)

where $\delta_{\mu,\nu}$ is the Kronecker delta function and
$\lambda^{2}=\textrm{max}\left(a_{c},r_{c}\times\textrm{max}(\hat{\sigma}_{\mu}^{2})\right)$
is a parameter of the regularization. We found that for Ansatz I, this method
produced the most stable and accurate results, with $a_{c}=r_{c}=10^{-4}$. On
the other hand, we found that for Ansatz II, discarding the smaller
eigenvalues with respect to a cutoff $\delta_{c}=10^{-9}$ and solving the
reduced linear system is a more suitable approach. For the simulations with
IBM’s Qasm simulator, the aforementioned regularization methods could not be
applied, as the positivity of the $M$ matrix is no longer preserved, with some
of its eigenvalues taking negative values. This was accounted for with a
variation of the diagonal shift method, which consists in adding a small
positive constant $\lambda$ to the diagonal terms of $M$ and expanding the
solution of Eq.(7) with respect to $\lambda$ as follows:

$\displaystyle\dot{\beta}$
$\displaystyle=M^{-1}V=(M+\lambda\mathbb{I}-\lambda\mathbb{I})^{-1}V$ (57)
$\displaystyle=(1-\lambda(M+\lambda\mathbb{I})^{-1})(M+\lambda\mathbb{I})^{-1}V$
$\displaystyle=\sum_{i=0}\lambda^{i}(M+\lambda\mathbb{I})^{-i}(M+\lambda\mathbb{I})^{-1}V$
$\displaystyle=(M+\lambda\mathbb{I})^{-1}V+\lambda(M+\lambda\mathbb{I})^{-1}(M+\lambda\mathbb{I})^{-1}V$
$\displaystyle+\lambda^{2}(M+\lambda\mathbb{I})^{-1}(M+\lambda\mathbb{I})^{-1}(M+\lambda\mathbb{I})^{-1}V+O(\lambda^{3})\,,$

or equivalently,

$\displaystyle\dot{\beta}$
$\displaystyle=\dot{\beta_{0}}+\lambda\dot{\beta_{1}}+\lambda^{2}\dot{\beta_{2}}+\lambda^{3}\dot{\beta_{3}}+...$
(58) $\displaystyle=\sum_{i=0}\lambda^{i}\dot{\beta_{i}}\,,$

with

$(M+\lambda\mathbb{I})\dot{\beta}_{0}=V\,,$ (59)

and

$(M+\lambda\mathbb{I})\dot{\beta}_{i}=\dot{\beta}_{i-1}\ ,i>0\,.$ (60)

For the LRQTE of a 2-site dissipative TFIM (Fig.6), we truncated the expansion
to second order in $\lambda$, and for the LRQTE of a 3-site dissipative TFIM
(Fig.7), we truncated to third order in $\lambda$.