A quantum algorithm for the direct estimation of the steady state of open quantum systems

Nathan Ramusat
Institute of Physics, Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland
Vincenzo Savona
Institute of Physics, Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland

Simulating the dynamics and the non-equilibrium steady state of an open quantum system are hard computational tasks on conventional computers. For the simulation of the time evolution, several efficient quantum algorithms have recently been developed. However, computing the non-equilibrium steady state as the long-time limit of the system dynamics is often not a viable solution, because of exceedingly long transient features or strong quantum correlations in the dynamics. Here, we develop an efficient quantum algorithm for the direct  estimation of averaged expectation values of observables on the non-equilibrium steady state,  thus bypassing the time integration of the master equation. The algorithm  encodes the vectorized representation of the density matrix on a quantum register, and makes use of quantum phase estimation to approximate the eigenvector associated to the zero eigenvalue of the generator of the system dynamics. We show that the output state of the algorithm allows to estimate expectation values of observables on the steady state. Away from critical points, where the Liouvillian gap scales as a power law of the system size, the quantum algorithm performs with exponential advantage compared to exact diagonalization.

§ INTRODUCTION

Open quantum systems are rapidly emerging as a major research field [Carusotto and Ciuti, 2013, Noh and Angelakis, 2017, Hartmann, 2016]. On the fundamental side, they can display new classes of universal physical properties such as dissipative phase transitions and topological phases, while applications can span from novel paradigms of quantum simulators [Rota et al., 2019, Rota and Savona, 2019] to the accurate modeling of noise in modern quantum computing platforms [Michael et al., 2016]. The numerical simulation of the dynamics of a many-body open quantum system is generally a hard computational task, as it combines the unitary evolution generated by the many-body Hamiltonian with the non-unitary evolution induced by the interaction with the environment. Similarly, computing the non-equilibrium steady state (NESS) – i.e. the state reached asymptotically in the long-time limit – is the computational analog of simulating the ground state of a many-body closed quantum system, and generally embodies an analogous computational challenge.

Current approaches [Weimer et al., 2019] to the solution of the Lindblad-Von-Neumann master equation for the density matrix [Breuer and Petruccione, 2002] – describing the case of a Markovian environment – include in addition to the exact numerical solution of the resulting differential equations, tensor-network methods [Werner et al., 2016, Verstraete et al., 2004, Zwolak and Vidal, 2004], projector [Nagy and Savona, 2018] and time-dependent variational [Nagy and Savona, 2019, Hartmann and Carleo, 2019, Vicentini et al., 2019, Yoshioka and Hamazaki, 2019] quantum Monte Carlo, various quantum trajectory approaches, as well as methods introducing various levels of approximation [Biella et al., 2018, Casteels et al., 2018, Verstraelen et al., 2020, Verstraelen and Wouters, 2018, Lange et al., 2018]. When the goal of the simulation is to know the NESS, integrating the system dynamics over long times is generally not an optimal strategy. In presence of critical slowing down, the time required to actually reach the NESS may be prohibitively large, and the transient dynamics may explore highly quantum correlated states – generally more difficult to simulate – even when the quantum correlations in the NESS are moderate. More specific approaches aiming directly at the simulation of the NESS without simulating the underlying dynamics have been developed, either in terms of variational tensor network [Kshetrimayum et al., 2017, Mascarenhas et al., 2015, Cui et al., 2015], and specific real-space decimation schemes [Finazzi et al., 2015]. The runtime of most of these approaches scales exponentially with the number of degrees of freedom, while approximate methods display power-law scaling but at the cost of a limited predictive power when in presence of significant quantum correlations in the system.

Quantum algorithms for the simulation of open quantum systems by computing the system dynamics have recently been proposed [Chenu et al., 2017, Wang et al., 2011, Barreiro et al., 2011, Kliesch et al., 2011, García-Pérez et al., 2020, Su and Li, 2020, Hu et al., 2020, Sweke et al., 2014, Sweke et al., 2015, Di Candia et al., 2015, Sweke et al., 2016, Childs and Li, 2017, Cleve and Wang, 2016, Liu et al., 2020]. The time-evolution of the Lindblad-Von-Neumann master equation is non-unitary, hence not suited for a direct implementation as a digital quantum simulation. While some special cases of system-environment interaction have been shown to translate into a unitary stochastic quantum evolution [Chenu et al., 2017], more generally the non-unitary dynamics is mapped onto a unitary dynamics on an appropriately dilated Hilbert space, which models an effective environment. Similarly, a time evolution quantum algorithm was proposed to compute the thermal equilibrium state of a system in interaction with a thermal bath [Terhal and DiVincenzo, 2000]. Recently, a hybrid quantum algorithm has been proposed for the estimation of the NESS, based on a variational ansatz for the density matrix [Yoshioka et al., 2020].

Here, we propose a different approach consisting in the direct estimation of the averaged expectation values of observables on the NESS, without requiring the integration of the system dynamics. This is achieved by representing the density matrix in vectorized form and mapping the steady-state condition onto a linear system of equations. The quantum algorithm then uses quantum phase estimation (QPE) and the HHL algorithm [Harrow et al., 2009] to solve the linear system. The output state is an estimate of the elements of the density matrix of the NESS, encoded in a vectorized representation onto the output state of a quantum register. The initial state is prepared leveraging the known spectral properties of the generator of the open-system dynamics, so to obtain a large overlap with the target output state. In this way, the success probability of a single QPE run is $O(1)$. It is shown that in cases in which the Liouvillian gap – i.e. the asymptotic relaxation rate towards the NESS – decreases as a power law of the system size, the algorithm runs with exponential speedup compared to exact diagonalization of the master equation on classical hardware. By using established methods for the quantum measurement of expectation values [Knill et al., 2007], is it finally shown how to directly estimate averages of the expectation value of physical observables on the NESS statistical ensemble.

The article is organized as follows. In Section 2 the master equation formalism is introduced and cast in a form suitable to be encoded on a quantum computer. In Section 3 the quantum algorithm is presented and an analysis of the success probability, gate cost and errors is performed. In Section 4 the formalism to estimate expectation values on the NESS is derived, and its accuracy characterized. Section 5 contains a numerical analysis of the algorithm applied to a simple driven-dissipative problem. The conclusions are presented in Section 6. Appendix A develops the oracular part of the algorithm for the specific case of a dissipative transverse Ising spin model.

§ FORMALISM

The dynamics of an open quantum system interacting with a Markovian environment is described by the Lindblad-Von Neumann master equation $\dot{\hat\rho}={\cal L}(\hat\rho)$ [Breuer and Petruccione, 2002], where $\hat\rho$ is the density operator, and the Liouvillian super-operator is expressed as (setting $\hbar=1$)

\begin{equation}
{\cal L}(\hat\rho)=-i\left[\hat H,\hat\rho\right]-\frac{1}{2}\sum_j\left(\left\{\hat A_j^\dagger\hat A_j,\hat\rho\right\}-2\hat A_j\hat\rho\hat A_j^\dagger\right)\,.
\label{Lindblad}
\end{equation}

Here, $\hat H$ is the system Hamiltonian and $\hat A_j$ are jump operators describing the transitions induced on the system by the environment. For the scope of this work, we will assume that both $\hat H$ and the $\hat A_j$'s are quasi-local operators, i.e. they are expressed as sums of tensor products of at most few local degrees of freedom of the system, so that they can be efficiently implemented in a quantum circuit [Lloyd, 1996]. The case of global jump operators, describing transitions between eigenstates of $\hat H$, can also be addressed, provided an efficient way of approximating these operators by a sequence of quantum gates exists. We assume that a unique non-equilibrium steady state (NESS) exists and is therefore characterized by the condition ${\cal L}(\hat\rho_{ss})=0$. The existence and uniqueness of the NESS can be established within rather general assumptions, and is in particular verified for problems defined in a finite-dimensional Hilbert space [Nigro, 2019], as it will always be the case on a quantum circuit implementation of the problem.

In order to encode the density matrix onto the state of a quantum register, we adopt the vectorized representation, whereby operators on the Hilbert space ${\cal H}$ are mapped onto vectors in the tensor product space ${\cal H}\otimes{\cal H}$. Given an orthonormal basis, a density operator $\hat\rho=\sum_{jk}\rho_{jk}|j\rangle\langle k|$ maps onto the vector $|\rho\rangle=\sum_{jk}\rho_{jk}|j\rangle\otimes|k\rangle$ (appropriately normalized). This essentially corresponds to building a vector whose components are the concatenated columns of the density matrix. Under this mapping, the Liouvillian becomes a linear operator acting on the space ${\cal H}\otimes{\cal H}$, defined as

\begin{align}
{\cal L}&=-i\left(I\otimes\hat H-\hat H^T\otimes I\right)\nonumber\\
&-\frac{1}{2}\sum_j\left(I\otimes\hat A_j^\dagger\hat A_j+\hat A_j^T\hat A_j^*\otimes I-2\hat A_j^*\otimes\hat A_j\right)\,,
\label{LindbladVec}
\end{align}

where $I$ denotes the identity matrix. Here, we will assume that the Hilbert space ${\cal H}$ of the problem has dimension $2^N$. Then the vector $|\rho\rangle$ is defined in a $2^{2N}$-dimensional space and can therefore be encoded onto a $2N$-qubit register.

The Liouvillian operator for a dissipative system is in general not Hermitian. Its eigenvalues $\lambda_j$ are in general complex, and the nonzero eigenvalues are characterized by $\mathrm{Re}(\lambda_j)<0$. The right and left eigenvectors associated to the zero-eigenvalue (here assumed unique) are respectively the NESS $\hat\rho_{ss}$ and the identity matrix, i.e. ${\cal L}|\rho_{ss}\rangle=0$ and $\langle I|{\cal L}=0$ [Breuer and Petruccione, 2002, Minganti et al., 2018]. The latter relation implies also ${\cal L}^\dagger|I\rangle=0$. We define a Hermitian operator in a $2^{2N+1}$-dimensional space as

\begin{equation}
M = \left(\begin{array}{ll}
0& {\cal L}\\
{\cal L}^\dagger &0
\end{array}\right)\,.
\end{equation}

From the spectral properties of the Liouvillian, we gather that the operator $M$ has two eigenvectors with zero eigenvalue, namely $|\eta_0\rangle=|0\rangle|I\rangle$ and $|\eta_1\rangle=|1\rangle|\rho_{ss}\rangle$. From the assumption that the NESS is unique, we conclude that these will be the only eigenvectors with zero eigenvalue of $M$. More generally, we denote with $|\eta_j\rangle$ and $\varphi_j$ the eigenvectors and the corresponding eigenvalues of $M$. Given the structure of $M$, for each eigenvalue $\varphi_j$ associated to the eigenvector $|\eta_j\rangle=|0\rangle|a_j\rangle+|1\rangle|b_j\rangle$, the vector $|\tilde\eta_j\rangle=|0\rangle|a_j\rangle-|1\rangle|b_j\rangle$ is also eigenvector of $M$ with eigenvalue $-\varphi_j$. The eigenvector condition also implies the relations ${\cal L}^\dagger{\cal L}|b_j\rangle=\varphi_j^2|b_j\rangle$ and ${\cal L}{\cal L}^\dagger|a_j\rangle=\varphi_j^2|a_j\rangle$, thus proving that the eigenvalues of $M$ coincide with those of the operators ${\cal L}^\dagger{\cal L}$ and ${\cal L}{\cal L}^\dagger$. Hence, the values $\varphi_j$ are the singular values of the Liouvillian ${\cal L}$ and, by applying Weyl's inequality for singular values, one can finally show that for the Liouvillian gap $g$ – i.e. the minimal value of $|\mathrm{Re}(\lambda_j)|$ among all nonzero eigenvalues – the inequality $g\le\min_{\varphi_j\ne0}|\varphi_j|$ holds. We conclude that the Liouvillian gap $g$ sets a lower bound to the spectral gap of the operator $M$. We will use this property in what follows, when discussing the computational complexity of the quantum algorithm.

To encode the operator $M$, the Liouvillian operator is expressed in terms of a Hermitian and an anti-Hermitian term as ${\cal L}={\cal L}_H-i{\cal L}_A$, with ${\cal L}_H^\dagger={\cal L}_H$ and ${\cal L}_A^\dagger={\cal L}_A$. The operator $M$ is then expressed as

\begin{equation}
M = X\otimes{\cal L}_H+Y\otimes{\cal L}_A\,,
\label{M}
\end{equation}

where $X$ and $Y$ are Pauli operators. If ${\cal L}_H$ and ${\cal L}_A$ can be efficiently encoded in terms of quasi-local operators, then $M$ can also be efficiently encoded. From Eq. (<ref>),

\begin{align}
{\cal L}_A&=\left(I\otimes\hat H-\hat H^T\otimes I\right)+\frac{i}{2}\sum_j\left(\hat A_j^*\otimes\hat A_j-\hat A_j^T\otimes\hat A_j^\dagger\right)\nonumber\\
{\cal L}_H&=\frac{1}{2}\sum_j\left(\hat A_j^*\otimes\hat A_j+\hat A_j^T\otimes\hat A_j^\dagger-I\otimes\hat A_j^\dagger\hat A_j-\hat A_j^T\hat A_j^*\otimes I\right)\,.
\label{LindbladAH}
\end{align}

§ QUANTUM ALGORITHM

The goal of the algorithm is to prepare an initial state having a large overlap with the degenerate zero eigenspace of $M$. QPE is then applied to estimate the projection of this initial state onto the zero eigenspace. The quantum circuit is illustrated in Fig. <ref>. It comprises an input-state preparation stage, a QPE stage, and a measurement stage.

@C=1em @R=0.2em
[-7.5mm]16     5H^⊗t                     ⋯        11         5FT^†     [-9mm]16  [7.5mm]16         5

H^⊗t                     ⋯                             FT^†

H^⊗t                     ⋯                             FT^†

|0⟩         H^⊗t                     ⋯                             FT^†

H^⊗t             7         ⋯                        FT^†

H^⊗t     6                 ⋯                        FT^†

5𝒫         5U^2^0         5U^2^1         ⋯        5U^2^t-1

ℱ         U^2^0         U^2^1         ⋯        U^2^t-1

ℱ         U^2^0         U^2^1         ⋯        U^2^t-1

|0⟩     ℱ         U^2^0         U^2^1         ⋯        U^2^t-1

ℱ         U^2^0         U^2^1         ⋯        U^2^t-1

ℱ         U^2^0         U^2^1         ⋯        U^2^t-1
Quantum circuit of the NESS solver. The three dashed lines indicate the stages where the states $|\psi_1\rangle$, $|\psi_2\rangle$, and $|\psi_3\rangle$, Eqs. (<ref>-<ref>), occur. The Quantum circuit of the gate $\cal P$ for the preparation of the initial state $|\xi\rangle$ for the second register is detailed in Fig. <ref>

\Qcircuit @C=1em @R=0.2em {
& \multigate{5}{\mathcal{P}} & \qw & \\
& \ghost{\mathcal{F}} & \qw &   \\
& \ghost{\mathcal{F}} & \qw &   \\
& \ghost{\mathcal{F}} & \qw &   \\
& \ghost{\mathcal{F}} & \qw &   \\
& \ghost{\mathcal{F}} & \qw &
\Qcircuit @C=1em @R=0.0em @!R {
& & \gate{H} & \ctrl{2} & \qw & \ctrl{3} & \qw & \qw & \cdots && \ctrl{4} & \gate{X} & \qw \\
& & & & & & & & & & & & \\
& & \qw & \gate{H} & \ctrl{4} & \qw & \qw & \qw & \cdots && \qw & \qw & \qw \\
& & \qw & \qw & \qw & \gate{H} & \ctrl{4} & \qw & \cdots && \qw & \qw & \qw &\rstick{N}\\
& & \qw & \qw & \qw & \qw & \qw & \qw & \cdots && \gate{H} & \ctrl{4} & \qw & & \\
& & & & & & & & & & & & \\
& & \qw & \qw & \targ & \qw & \qw & \qw & \cdots && \qw & \qw & \qw \\
& & \qw & \qw & \qw & \qw & \targ & \qw & \cdots && \qw & \qw & \qw &\rstick{N}\\
& & \qw & \qw & \qw & \qw & \qw & \qw & \cdots && \qw & \targ & \qw
\gategroup{7}{13}{9}{13}{.7em}{\}}
\gategroup{3}{13}{5}{13}{.7em}{\}}

Quantum circuit $\cal P$ for the preparation of the initial state $|\xi\rangle$ for the second register. The circuit takes as input the $(2N+1)$-qubit state $|\mathbf{0}\rangle$ and outputs the state $|\xi\rangle$ in Eq. (<ref>). Controlled-Hadamard gates are implemented following the general procedure for arbitrary controlled gates [Nielsen and Chuang, 2016, Barenco et al., 1995], and controlled-arbitrary gates are available within current superconducting circuit technology [Lacroix et al., 2020].

The QPE algorithm uses two registers. The first one is a $t$-qubit register that will encode a $t$-bit integer estimate of the eigenvalue. The second one encodes the corresponding eigenvector and therefore contains, in our formulation, $2N+1$ qubits. In its standard formulation [Nielsen and Chuang, 2016], QPE takes as an input the state $|\mathbf{0}\rangle|\eta_j\rangle$, where $|\eta_j\rangle$ is an eigenstate of a unitary operator $U$ with eigenvalue $e^{2\pi i\varphi_j}$, assuming $\varphi_j\in[0,1]$. If $2^t\varphi_j$ is a $t$-bit integer, then QPE outputs the state $|\varphi_j\rangle|\eta_j\rangle$. If $2^t\varphi_j$ is not an integer, then QPE will produce with high probability the state $|\tilde\varphi_j\rangle|\eta_j\rangle$, where $2^t\tilde\varphi_j$ is the best $t$-bit integer estimate of $2^t\varphi_j$ from below.

For the present case, we set $U=e^{2\pi iMt_0}$, where the real parameter $t_0$ is set so that the whole spectrum of $Mt_0$ is included in the $[0,1]$ interval. The algorithm assumes that an oracle circuit is available to efficiently perform controlled-$U^{2^j}$ operations.

Differently from the standard QPE scheme, here the eigenvalue is known and the algorithm is used to compute the corresponding eigenvector. To this purpose, the input state $|\mathbf{0}\rangle|\xi\rangle$ of the QPE stage of the algorithm must be chosen in such a way that $|\xi\rangle$ has a significant overlap with the zero-eigenspace of $M$ spanned by $|\eta_0\rangle$ and $|\eta_1\rangle$. Here we set

\begin{equation}
\label{xi}
\end{equation}

This choice is justified by the procedure to estimate the expectation values of observables that will be introduced below. Insight in the structure of the state $|I\rangle$ shows immediately that it can be expressed as a product of Bell states, where each Bell state involves one qubit from the row index and one from the column index of the density matrix. A circuit on the second register, to prepare the input state (<ref>) when input with the initial state $|\mathbf{0}\rangle$, is shown in Fig. <ref>. Controlled Hadamard gates can be implemented with a limited number of gates following the prescriptions for arbitrary controlled operations [Nielsen and Chuang, 2016, Barenco et al., 1995], and controlled-arbitrary gates are available within current superconducting circuit technology [Lacroix et al., 2020].

Expanding $|\xi\rangle$ on the basis of the eigenvectors $|\eta_j\rangle$ of $M$, the state of the full circuit following the preparation stage is

\begin{align}
|\psi_1\rangle &= |\mathbf{0}\rangle|\xi\rangle \nonumber\\
&=  |\mathbf{0}\rangle\frac{|0\rangle|I\rangle+|1\rangle|\mathbf{0}\rangle}{\sqrt{2}}\nonumber\\
&=  \frac{|\mathbf{0}\rangle|\eta_0\rangle}{\sqrt{2}}+\frac{1}{\sqrt{2}}\sum_{j\ne0}c_{j}|\mathbf{0}\rangle|\eta_j\rangle\nonumber\\
&=  \frac{|\mathbf{0}\rangle|\eta_0\rangle+c_1|\mathbf{0}\rangle|\eta_1\rangle}{\sqrt{2}}+\frac{1}{\sqrt{2}}\sum_{j\ne0,1}c_j|\mathbf{0}\rangle|\eta_j\rangle\,,
\label{psi1}
\end{align}

where we have indicated explicitly the expansion of the state $|1\rangle|\mathbf{0}\rangle$ on the basis of eigenstates $|\eta_j\rangle$ of $M$, and accounted for the fact that states $|1\rangle|\mathbf{0}\rangle$ and $|\eta_0\rangle=|0\rangle|I\rangle$ are orthogonal. The first term of the sum in the last line is the one that will be untouched by the QPE, as both $|\eta_0\rangle$ and $|\eta_1\rangle$ are eigenvectors of $M$ associated to the zero eigenvalue.

In general, the other eigenvalues of $M$ will be such that $2^t\varphi_j$ is not a $t$-bit integer. Then, the QPE stage of the circuit applied to the state $|\mathbf{0}\rangle|\eta_j\rangle$ will output the state $\sum_{k=0}^{2^t-1}\alpha_k^{(j)}|\mathbf{k}\rangle|\eta_j\rangle$, so that the state of the circuit at the output of the QPE stage is

\begin{equation}
|\psi_2\rangle = \frac{|\mathbf{0}\rangle|\eta_0\rangle+c_1|\mathbf{0}\rangle|\eta_1\rangle}{\sqrt{2}}+\frac{1}{\sqrt{2}}\sum_{j\ne0,1}\sum_{k=0}^{2^t-1}c_j\alpha_k^{(j)}|\mathbf{k}\rangle|\eta_j\rangle\,,
\label{psi2}
\end{equation}

\begin{equation}
\alpha_k^{(j)}=\frac{1}{2^t}\frac{1-e^{2\pi i(2^t\varphi_j-k)}}{1-e^{2\pi i(\varphi_j-k/2^t)}}\,.
\end{equation}

After the QPE stage, the first register is measured in the computational basis, and the algorithm is successful if the eigenvalue $\mathbf{0}$ is measured. Then, the output state is (up to a normalization constant of order 1)

\begin{equation}
|\psi_3\rangle = |\mathbf{0}\rangle\left[\frac{|\eta_0\rangle+c_1|\eta_1\rangle}{\sqrt{2}}+\frac{1}{\sqrt{2}}\sum_{j\ne0,1}c_j\alpha_0^{(j)}|\eta_j\rangle\right]\,,
\label{psi3}
\end{equation}

where $c_1\in\mathbb{R}$ because it is the projected amplitude of the density matrix onto its first diagonal matrix element.
The first term in the sum between brackets in expression (<ref>) is the result being sought, while the remaining term is an unwanted error proportional to the eigenstates $|\eta_j\rangle$ of $M$ other than those with zero eigenvalues.

§.§ Success probability and runtime

From (<ref>), the probability of measuring zero on the first register is

\begin{align}
\label{p0}
\end{align}

where $p_e=1/2\sum_{j\ne0,1}|c_j|^2|\alpha_0^{(j)}|^2$ is the probability of an error in the QPE, whereby the zero eigenvalue is measured but the second register is projected onto one of the states $|\eta_j\rangle$ with $j\ne0,1$.
As $p_0>1/2$, the number of executions required to successfully measure the zero eigenvalue is $O(1)$ and does not scale with $N$ or $t$.

In order to estimate the contribution of the spurious term in (<ref>), it is useful to study the behavior of $|\alpha_0^{(j)}|^2$. Simple algebra leads to

\begin{equation}
\label{sinsin}
\end{equation}

The function (<ref>) is plotted in Fig. <ref>(a). As already pointed out, the value of $t_0$ is set in order to have $\varphi_j<1$. The value of $|\alpha_0^{(j)}|^2$ can be large only when $\varphi_j$ approaches an integer value, and the value of $t_0$ can be set such that the only relevant case is when $\varphi_j\sim0$. For $\varphi_j\sim0$ we can expand the denominator and obtain (setting $x=\pi2^t\varphi_j$) $|\alpha_0^{(j)}|^2\simeq\sinc^2(x)$, the squared sine cardinal function. We can use the property $\sinc^2(x)\le x^{-2}$ to set a bound on the error probability

\begin{equation*}
\end{equation*}

where $g$ is the Liouvillian spectral gap and we have used the inequality previously established. Then, in order to have a bound probability of error $p_e<\epsilon_p$, one must set

\begin{equation*}
t>\left\lceil\log\left(\frac{1}{\sqrt{2}\pi g}\right)+\log\left(\frac{1}{\epsilon_p}\right)\right\rceil\,.
\label{nqubit}
\end{equation*}

(a) The probability $|\alpha_0^{(j)}|^2$ plotted as a function of $\varphi_j$ for different values of $t$. (b) The accuracy $1-{\cal F}$ of the quantum algorithm as a function of $t$, as computed for the simple model of one spin in presence of an external field $h$ along the $x$-direction and relaxation along the $z$-direction.

To estimate the gate cost of the algorithm and its scaling with the system size $N$, we first observe that the initial state preparation ${\cal P}$ is made of $2N+2$ gates. The Hadamard stage $H^{\otimes t}$ contributes for $t$ gates, while the QFT stage is implemented with $O(t^2)$ gates if the textbook implementation is used, while approximate versions of the QFT can be implemented with $O(t\log(t))$ gates [Nam et al., 2020, Cleve and Watrous, 2000]. The oracle part of QPE amounts to the application of a controlled $U^{2^t}=e^{2\pi i 2^tMt_0}$ to the $(2N+1)$-qubit register. The control is implemented at no additional gate cost if a Suzuki-Trotter expansion [Childs et al., 2021] is used, as shown explicitly in Appendix A for the case of the dissipative transverse-field Ising model.

A general physical setting is that of $k$-local Hamiltonians $\hat H=\sum_{\gamma}\alpha_{\gamma}\hat H_{\gamma}$ where each $\hat H_{\gamma}$ acts in a non trivial way on $k\sim O(1)$ spins only. We additionally assume that the jump operators $\hat A_j$ in Eq. (<ref>) are also quasi-local operators. Then, the vectorized Liouvillian Eq. (<ref>) is $k$-local and, from Eq. (<ref>), it follows that the effective Hamiltonian $M$ is $(k+1)$-local on an effective system of size $2N+1$. An example is provided by the $3$-local expression of $M$ in the case of the transverse Ising model with local spin relaxation (<ref>) derived in Appendix A. Optimal strategies for the approximate quantum simulation of Hamiltonian evolution over a time $\tau$ to error $\epsilon_\tau$ have been developed, adopting both Suzuki-Trotter expansion or alternative approaches [Childs et al., 2021]. The Suzuki-Trotter approach has the advantage of not requiring ancilla qubits. Assuming a $p_{th}$-order Suzuki-Trotter expansion and the most general $k$-local operator $M$, the best gate cost for the implementation of a $C-e^{i2\pi t_0M\tau}$ gate is estimated [Childs et al., 2021] as $O((2N+1)^k\tau^{(1+1/p)}/\epsilon_\tau^{1/p})$ up to the 1-norm of the operator $2\pi t_0||M||$ which in the present case is assumed to be $O(1)$. As noted above, here $\tau=2^t$, which leads us to an estimation of the gate cost $O((2N+1)^k/( g^{1+1/p}\epsilon_\tau^{1/p}\epsilon_p^{1/2})$ for the simulation of the vectorized NESS of a $N$-spin system characterized by a Liouvillian gap $g$, with probability of error $p_e<\epsilon_p$.
The $(2N+1)^k$-dependence is determined by the number of gates required for the most general $k$-local interaction term in the Hamiltonian. However, for models with limited range interactions, the corresponding factor typically depends linearly on $N$ [Tacchino et al., 2020], as for example in the case of the transverse Ising model discussed in Appendix A. The overall gate cost of the present algorithm is therefore

\begin{equation}
\end{equation}

where $\epsilon_\tau$ and $\epsilon_p$ are respectively the Suzuki-Trotter error and the QPE success error.

Choosing a sufficiently high order $p$, the gate cost approaches $O((2N+1)^k/(g\epsilon_p^{1/2}))$, i.e. linear in the inverse Liouvillian gap and independent of the Suzuki-Trotter error $\epsilon_\tau$. The inverse gap for most model Hamiltonians and couplings to the environment scales as a power law of the system size, except at critical points where it can scale exponentially [Žnidarič, 2015]. The required number of qubits $t$ in the first register, according to Eq. (<ref>), is $t\sim O(\log(N))$ in the case of power-law Liouvillian gap, and $t\sim O(N)$ otherwise. Assuming a power-law dependence $1/g=\mathrm{poly}(N)$, the overall gate cost of the algorithm scales with the system size as $O(\mathrm{poly}(N))$, thus achieving a quantum advantage with respect to an exact calculation with a classical algorithm. The scaling with $N$ is instead not favorable whenever the inverse gap scales exponentially with system size. This can happen in some systems when approaching the critical point of a dissipative phase transition, in a way analogous to the simulation of the ground state of a Hamiltonian close to a quantum phase transition using either QPE [Nielsen and Chuang, 2016] or a variational quantum eigensolver [McClean et al., 2016, Peruzzo et al., 2014].

Previous approaches to the quantum simulation of open quantum systems [Chenu et al., 2017, Wang et al., 2011, Barreiro et al., 2011, Kliesch et al., 2011, García-Pérez et al., 2020, Su and Li, 2020, Hu et al., 2020, Sweke et al., 2014, Sweke et al., 2015, Di Candia et al., 2015, Sweke et al., 2016, Childs and Li, 2017, Cleve and Wang, 2016, Liu et al., 2020] have focused on algorithms to simulate the time evolution of the system. The scope of the present work is different, in that we propose an algorithm to solve directly the algebraic equation ${\cal L}(\hat\rho_{ss})=0$ that defines the density matrix associated to the NESS. Computing the NESS through time evolution would require a time $T$ long enough to asymptotically approximate the NESS. This time scale is naturally defined as the inverse Liouvillian gap $T\sim1/g$. State-of-the-art approaches [Cleve and Wang, 2016] achieve a gate cost $O(T\mathrm{polylog}(T/\epsilon)$, for a required accuracy $\epsilon$. This is comparable to our approach which scales as $O(T^{1+1/p}/\epsilon^{1/p})$ for a $p_{th}$-order Suzuki-Trotter expansion and approaches a complexity linear in $T$ in the limit of large $p$. The merit of the present approach is that it bypasses the computation of the real time evolution altogether. The transient to the NESS is often characterized by short time features due to the Hamiltonian part of the evolution, that would require a small integration time step, possibly leading to a system dependent computational overhead. Computing directly the density matrix of the NESS avoids this issue and is therefore the method to be preferred if only averaged expectation values of observables on the steady state are needed. A further distinctive feature of the present approach is the fact that its output state encodes directly the elements of the density matrix characterizing the steady state. Most previous approaches instead [Wang et al., 2011, Barreiro et al., 2011, Kliesch et al., 2011, García-Pérez et al., 2020, Su and Li, 2020, Hu et al., 2020, Sweke et al., 2014, Sweke et al., 2015, Sweke et al., 2016, Childs and Li, 2017, Cleve and Wang, 2016] simulate the quantum channel associated to the open system being studied, typically through appropriate dilation schemes. To estimate the ensemble averaged expectation values of observables, these approaches thus additionally require a possibly large number of runs in order to sample over the statistical ensemble of pure states characterizing the quantum channel. If the goal of the quantum simulation is to estimate ensemble-averaged expectation values, then the present algorithm will likely bring a significant computational advantage. A similar goal is pursued in Ref. [Di Candia et al., 2015] where the time evolution of correlators is simulated, although the perturbative scheme used for the dissipative dynamics results in a computational cost scaling exponentially with time, thus inappropriate to extrapolate the NESS in the long-time limit.

§ ESTIMATE OF EXPECTATION VALUES

Once the first register has been measured in the state $|\mathbf{0}\rangle$, it is possible to estimate quantum mechanical expectation values from the density matrix $\hat\rho_{ss}$. For a given observable $\hat O$, the expectation value is computed as $\langle\hat O\rangle=\mbox{Tr}(\hat O\hat\rho_{ss})/\mbox{Tr}(\hat\rho_{ss})$. As before, we assume that the Hermitian operator $\hat O$ can be efficiently encoded in terms of one- and two-qubit operations, which is generally true for few-point correlations of quasi-local observables. In vectorized form, we have for a general density matrix $|(\hat O\hat\rho)\rangle=I\otimes\hat O|\rho\rangle$, and the super-operator ${\cal O}=I\otimes\hat O$ is also Hermitian. The trace average is then expressed as a matrix element

\begin{equation}
\mbox{Tr}(\hat O\hat\rho_{ss})=\langle I|{\cal O}|\rho_{ss}\rangle=\langle\rho_{ss}|{\cal O}|I\rangle\,.
\label{trace}
\end{equation}

In order to estimate an averaged expectation value (<ref>) from the output state of the algorithm, we define the Hermitian operator

\begin{equation}
{\cal Q} = X\otimes {\cal O}= \left(\begin{array}{ll}
0& {\cal O}\\
{\cal O} &0
\end{array}\right)\,,
\end{equation}

acting on the $(2N+1)$-qubit second register. Then,

\begin{equation}
\mathrm{Tr}(\hat O\hat\rho_{ss})=\langle\eta_0|{\cal Q}|\eta_1\rangle=\langle\eta_1|{\cal Q}|\eta_0\rangle\,,
\label{trace2}
\end{equation}

while, due to the structure of the operator $\cal Q$, we have $\langle\eta_0|{\cal Q}|\eta_0\rangle=\langle\eta_1|{\cal Q}|\eta_1\rangle=0$.

The trace average can then be estimated as the expectation value of $I\otimes{\cal Q}$ on the output state $|\psi_3\rangle$

\begin{align}
\langle\psi_3|I\otimes{\cal Q}|\psi_3\rangle&=c_1\mathrm{Tr}(\hat O\hat\rho_{ss})\label{exp}\\
&+\sum_{j\ne0,1}\mathrm{Re}\left[c_j\alpha_0^{(j)}\left(\langle\eta_0|{\cal Q}|\eta_j\rangle+c_1\langle\eta_1|{\cal Q}|\eta_j\rangle\right)\right]+\frac{1}{2}\sum_{j,l\ne0,1}c^*_jc_l\alpha_0^{(j)*}\alpha_0^{(l)}\langle\eta_j|{\cal Q}|\eta_l\rangle\nonumber\,,
\end{align}

where we used the facts that the expectation value of ${\cal O}$ is real-valued, $\langle\eta_0|{\cal Q}|\eta_1\rangle=\langle\eta_1|{\cal Q}|\eta_0\rangle$, and $c_1\in\mathbb{R}$. In particular, by setting $\hat O=I$, the measurement returns the quantity $c_1\mbox{Tr}(\hat\rho_{ss})$, which can be then used to eliminate the factor $c_1$ and obtain the expectation value $\langle\hat O\rangle=\mbox{Tr}(\hat O\hat\rho_{ss})/\mbox{Tr}(\hat\rho_{ss})$. The second and third term on the right hand side of Eq. (<ref>) are the spurious terms resulting from the finite number of qubits $t$ in the QPE. They are respectively of order $O(|\alpha_0^{(j)}|)$ and $O(|\alpha_0^{(j)}|^2)$. By setting $t_0$ as discussed in the previous Section, these terms therefore decrease exponentially in $t$, as also seen in the example presented below.

Eqs. (<ref>) and (<ref>) prove that the problem of computing the trace average of $\hat O$ maps onto that of estimating the expectation value of a Hermitian operator onto the output quantum state of the present algorithm. This estimate is achieved through methods for the quantum measurement of expectation values [Knill et al., 2007, Roggero and Baroni, 2020]. The accuracy of these methods scales with the number of measurements $N_m$ performed on the output state after, thus requiring $N_m$ successful runs of the algorithm. In order to achieve a given accuracy $\epsilon_m$, a direct operator averaging method requires $N_m=O(1/\epsilon_m^2)$ and no additional gates, while methods based on the QPE can achieve $N_m=O(1/\epsilon_m)$ but at an additional algorithmic cost [Knill et al., 2007, Roggero and Baroni, 2020].

An important question concerns the normalization of the output state encoding the vectorized density matrix, and how it affects the error in estimating expectation values. The appropriately normalized state is written in terms of the elements of the density matrix as $|\rho\rangle=(1/{\cal N})\sum_{jk}\rho_{jk}|j\rangle\otimes|k\rangle$, where ${\cal N}=\sqrt{\sum_{jk}|\rho_{jk}|^2}=\sqrt{P}$ is the Frobenius norm of the density matrix and $P=\mathrm{Tr}(\hat\rho^2)$ the purity of the NESS. If we assume the output state $|\psi_3\rangle$ to be normalized, then the expectation value in Eq. (<ref>) is actually $\langle\psi_3|I\otimes{\cal Q}|\psi_3\rangle=(c_1/P)\mathrm{Tr}(\hat O\hat\rho_{ss})$. In the case of a highly mixed NESS, the large $1/P$-factor amplifies the error in the estimation of $\langle\hat O\rangle$. For a quantum system described in a $2^N$-dimensional Hilbert space, the inequality $1\le1/P\le2^{N}$ holds. The lower bound corresponds to pure states, whereas the upper bound is set by the maximally mixed state corresponding to $\hat\rho=I/2^N$. In practice, physical models of interest usually have poorly entropic steady states, for which $1/P$ is either bound or scales sub-extensively close to a dissipative phase transition [Rota et al., 2017, Rota et al., 2019].

§ RESULTS

Estimated expectation values of (a) $\hat\sigma_y$ and (b) $\hat\sigma_z$, as computed from Eq. (<ref>). The horizontal lines correspond to the values obtained from the exact NESS. The insets show the relative errors $\delta=|\langle\hat\sigma_\alpha\rangle_\mathrm{est}-\langle\hat\sigma_\alpha\rangle|/|\langle\hat\sigma_\alpha\rangle|$ ($\alpha=y,z$) and the black line is a guide to the eye proportional to $2^{-t}$.

We present here a simple numerical test of the quantum algorithm by simulating the quantum circuit for an elementary model of an open quantum system. We consider a spin 1/2 subject to an external magnetic field $h$ along the $x$ direction and to a decay process along the z-axis. The system Hamiltonian is simply $\hat H=h\hat\sigma_x$ and the density matrix obeys the Lindblad-Von Neumann master equation

\begin{equation}
\dot{\hat\rho}=-i\left[\hat H,\hat\rho\right]-\frac{1}{2}\left(\hat\sigma^{+}\hat\sigma^{-}\hat\rho+\hat\rho\hat\sigma^{+}\hat\sigma^{-}-2\hat\sigma^{-}\hat\rho\hat\sigma^{+}\right)\,.
\end{equation}

The Liouvillian gap for this simple model is independent of $h$ and given by $g=1/2$. We therefore set $t_0=1/5$. The circuit in Figs. <ref> and  <ref> was simulated by explicitly computing the corresponding unitary operator, for a number of qubits in the first register up to $t=10$ and three qubits in the second register.

For varying values of $h$ and $t$, we computed the NESS density matrix $\hat\rho_{ss}^{(q)}$ resulting from the quantum algorithm, and compared it to the exact NESS $\hat\rho_{ss}$ by evaluating the fidelity

\begin{equation}
{\cal F}=\left(\mathrm{tr}\sqrt{\sqrt{\hat\rho_{ss}}~\hat\rho_{ss}^{(q)}\sqrt{\hat\rho_{ss}}}\right)^2\,.
\end{equation}

Results are displayed in Fig. <ref>(b). The accuracy of the quantum algorithm scales as $2^{-2t }$, as expected from Eqs. (<ref>) and (<ref>).

In Fig. <ref> the expectation values of $\hat\sigma_y$ and $\hat\sigma_z$, estimated according to Eq. (<ref>), are plotted as a function of the number of qubits $t$. The corresponding errors with respect to the exact expectation values are plotted in the insets. The leading term linear in $\alpha_0^{(j)}$ in Eq. (<ref>) results in a relative error scaling as $\delta\sim2^{-t}$, again proving the efficiency of the algorithm.

§ CONCLUSIONS

We have developed a quantum algorithm for a fault-tolerant quantum computer, that directly estimates the vectorized representation of the NESS density matrix of an open quantum system, without the need to integrate the master equation for the time evolution of the density operator. The algorithm leverages QPE to find the Liouvillian eigenvector associated to the null eigenvalue, thus assuming the existence of a unique NESS. Once this eigenvector is found, we show how to efficiently compute the averaged expectation value of an observable onto the NESS. The scaling of the present quantum algorithm is polynomial in the number of qubits $N$ encoding the quantum degrees of freedom of the system, thus providing quantum advantage when compared to exact diagonalization on classical computer hardware. The algorithm is designed to run on a fault-tolerant quantum computer. As such, we expect it to still execute correctly in a shallow circuit setting where the overall error rate is contained. The question about the scaling of its accuracy, when executed on near-term noisy quantum hardware, remains however open.

With the availability of scalable fault-tolerant quantum hardware, the present algorithm will enable modeling open quantum systems of unprecedented size, opening a new way in the study of the dissipation and decoherence, including among others, phenomena such as dissipative phase transitions or the influence of the environment on quantum hardware itself.

We acknowledge enlightening discussions with Fabrizio Minganti.

[Carusotto and Ciuti, 2013]
Iacopo Carusotto and Cristiano Ciuti.
Quantum fluids of light.
Reviews of Modern Physics, 850 (1):0 299–366,
feb 2013.
ISSN 0034-6861.
URL <https://link.aps.org/doi/10.1103/RevModPhys.85.299>.

[Noh and Angelakis, 2017]
Changsuk Noh and Dimitris G Angelakis.
Quantum simulations and many-body physics with light.
Reports on Progress in Physics, 800 (1):0
016401, jan 2017.
ISSN 0034-4885.

[Hartmann, 2016]
Michael J Hartmann.
Quantum simulation with interacting photons.
Journal of Optics, 180 (10):0 104005, sep
URL <https://doi.org/10.1088/2040-8978/18/10/104005>.

[Rota et al., 2019]
Riccardo Rota, Fabrizio Minganti, Cristiano Ciuti, and Vincenzo Savona.
Quantum critical regime in a quadratically driven nonlinear photonic
Phys. Rev. Lett., 122:0 110405, Mar 2019.
URL <https://link.aps.org/doi/10.1103/PhysRevLett.122.110405>.

[Rota and Savona, 2019]
Riccardo Rota and Vincenzo Savona.
Simulating frustrated antiferromagnets with quadratically driven qed
Phys. Rev. A, 100:0 013838, Jul 2019.
URL <https://link.aps.org/doi/10.1103/PhysRevA.100.013838>.

[Michael et al., 2016]
Marios H. Michael, Matti Silveri, R. T. Brierley, Victor V. Albert, Juha
Salmilehto, Liang Jiang, and S. M. Girvin.
New class of quantum error-correcting codes for a bosonic mode.
Phys. Rev. X, 6:0 031006, Jul 2016.
URL <https://link.aps.org/doi/10.1103/PhysRevX.6.031006>.

[Weimer et al., 2019]
Hendrik Weimer, Augustine Kshetrimayum, and Román Orús.
Simulation methods for open quantum many-body systems.
arXiv e-prints, art. arXiv:1907.07079, July 2019.
URL <https://arxiv.org/abs/1907.07079>.

[Breuer and Petruccione, 2002]
H. P. Breuer and F. Petruccione.
The theory of open quantum systems.
Oxford University Press, Great Clarendon Street, 2002.

[Werner et al., 2016]
A. H. Werner, D. Jaschke, P. Silvi, M. Kliesch, T. Calarco, J. Eisert, and
S. Montangero.
Positive tensor network approach for simulating open quantum
many-body systems.
Phys. Rev. Lett., 116:0 237201, Jun 2016.
URL <https://link.aps.org/doi/10.1103/PhysRevLett.116.237201>.

[Verstraete et al., 2004]
F. Verstraete, J. J. García-Ripoll, and J. I. Cirac.
Matrix product density operators: Simulation of finite-temperature
and dissipative systems.
Phys. Rev. Lett., 93:0 207204, Nov 2004.
URL <https://link.aps.org/doi/10.1103/PhysRevLett.93.207204>.

[Zwolak and Vidal, 2004]
Michael Zwolak and Guifré Vidal.
Mixed-state dynamics in one-dimensional quantum lattice systems: A
time-dependent superoperator renormalization algorithm.
Phys. Rev. Lett., 93:0 207205, Nov 2004.
URL <https://link.aps.org/doi/10.1103/PhysRevLett.93.207205>.

[Nagy and Savona, 2018]
Alexandra Nagy and Vincenzo Savona.
Driven-dissipative quantum monte carlo method for open quantum
Phys. Rev. A, 97:0 052129, May 2018.
URL <https://link.aps.org/doi/10.1103/PhysRevA.97.052129>.

[Nagy and Savona, 2019]
Alexandra Nagy and Vincenzo Savona.
Variational quantum monte carlo method with a neural-network ansatz
for open quantum systems.
Phys. Rev. Lett., 122:0 250501, Jun 2019.
URL <https://link.aps.org/doi/10.1103/PhysRevLett.122.250501>.

[Hartmann and Carleo, 2019]
Michael J. Hartmann and Giuseppe Carleo.
Neural-network approach to dissipative quantum many-body dynamics.
Phys. Rev. Lett., 122:0 250502, Jun 2019.
URL <https://link.aps.org/doi/10.1103/PhysRevLett.122.250502>.

[Vicentini et al., 2019]
Filippo Vicentini, Alberto Biella, Nicolas Regnault, and Cristiano Ciuti.
Variational neural-network ansatz for steady states in open quantum
Phys. Rev. Lett., 122:0 250503, Jun 2019.
URL <https://link.aps.org/doi/10.1103/PhysRevLett.122.250503>.

[Yoshioka and Hamazaki, 2019]
Nobuyuki Yoshioka and Ryusuke Hamazaki.
Constructing neural stationary states for open quantum many-body
Phys. Rev. B, 99:0 214306, Jun 2019.
URL <https://link.aps.org/doi/10.1103/PhysRevB.99.214306>.

[Biella et al., 2018]
Alberto Biella, Jiasen Jin, Oscar Viyuela, Cristiano Ciuti, Rosario Fazio, and
Davide Rossini.
Linked cluster expansions for open quantum systems on a lattice.
Phys. Rev. B, 97:0 035103, Jan 2018.
URL <https://link.aps.org/doi/10.1103/PhysRevB.97.035103>.

[Casteels et al., 2018]
Wim Casteels, Ryan M. Wilson, and Michiel Wouters.
Gutzwiller monte carlo approach for a critical dissipative spin
Phys. Rev. A, 97:0 062107, Jun 2018.
URL <https://link.aps.org/doi/10.1103/PhysRevA.97.062107>.

[Verstraelen et al., 2020]
Wouter Verstraelen, Riccardo Rota, Vincenzo Savona, and Michiel Wouters.
Gaussian trajectory approach to dissipative phase transitions: The
case of quadratically driven photonic lattices.
Phys. Rev. Research, 2:0 022037, May 2020.
URL <https://link.aps.org/doi/10.1103/PhysRevResearch.2.022037>.

[Verstraelen and Wouters, 2018]
Wouter Verstraelen and Michiel Wouters.
Gaussian quantum trajectories for the variational simulation of open
quantum-optical systems.
Applied Sciences, 80 (9), 2018.
ISSN 2076-3417.
URL <https://www.mdpi.com/2076-3417/8/9/1427>.

[Lange et al., 2018]
Florian Lange, Zala Lenarčči čč, and Achim Rosch.
Time-dependent generalized gibbs ensembles in open quantum systems.
Phys. Rev. B, 97:0 165138, Apr 2018.
URL <https://link.aps.org/doi/10.1103/PhysRevB.97.165138>.

[Kshetrimayum et al., 2017]
Augustine Kshetrimayum, Hendrik Weimer, and Román Orús.
A simple tensor network algorithm for two-dimensional steady states.
Nature Communications, 80 (1):0 1291, 2017.
ISSN 2041-1723.
URL <https://doi.org/10.1038/s41467-017-01511-6>.

[Mascarenhas et al., 2015]
Eduardo Mascarenhas, Hugo Flayac, and Vincenzo Savona.
Matrix-product-operator approach to the nonequilibrium steady state
of driven-dissipative quantum arrays.
Phys. Rev. A, 92:0 022116, Aug 2015.
URL <https://link.aps.org/doi/10.1103/PhysRevA.92.022116>.

[Cui et al., 2015]
Jian Cui, J. Ignacio Cirac, and Mari Carmen Bañuls.
Variational matrix product operators for the steady state of
dissipative quantum systems.
Phys. Rev. Lett., 114:0 220601, Jun 2015.
URL <https://link.aps.org/doi/10.1103/PhysRevLett.114.220601>.

[Finazzi et al., 2015]
S. Finazzi, A. Le Boité, F. Storme, A. Baksic, and C. Ciuti.
Corner-space renormalization method for driven-dissipative
two-dimensional correlated systems.
Phys. Rev. Lett., 115:0 080604, Aug 2015.
URL <https://link.aps.org/doi/10.1103/PhysRevLett.115.080604>.

[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:0 140403, Apr 2017.
URL <https://link.aps.org/doi/10.1103/PhysRevLett.118.140403>.

[Wang et al., 2011]
Hefeng Wang, S. Ashhab, and Franco Nori.
Quantum algorithm for simulating the dynamics of an open quantum
Phys. Rev. A, 83:0 062317, Jun 2011.
URL <https://link.aps.org/doi/10.1103/PhysRevA.83.062317>.

[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, 4700 (7335):0 486–491, Feb 2011.
ISSN 1476-4687.
URL <https://doi.org/10.1038/nature09801>.

[Kliesch et al., 2011]
M. Kliesch, T. Barthel, C. Gogolin, M. Kastoryano, and J. Eisert.
Dissipative quantum church-turing theorem.
Phys. Rev. Lett., 107:0 120501, Sep 2011.
URL <https://link.aps.org/doi/10.1103/PhysRevLett.107.120501>.

[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, 60 (1):0 1, 2020.
ISSN 2056-6387.
URL <https://doi.org/10.1038/s41534-019-0235-y>.

[Su and Li, 2020]
Hong-Yi Su and Ying Li.
Quantum algorithm for the simulation of open-system dynamics and
Phys. Rev. A, 101:0 012328, Jan 2020.
URL <https://link.aps.org/doi/10.1103/PhysRevA.101.012328>.

[Hu et al., 2020]
Zixuan Hu, Rongxin Xia, and Sabre Kais.
A quantum algorithm for evolving open quantum dynamics on quantum
computing devices.
Scientific Reports, 100 (1):0 3301, 2020.
ISSN 2045-2322.
URL <https://doi.org/10.1038/s41598-020-60321-x>.

[Sweke et al., 2014]
Ryan Sweke, Ilya Sinayskiy, and Francesco Petruccione.
Simulation of single-qubit open quantum systems.
Phys. Rev. A, 90:0 022331, Aug 2014.
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:0 062308, Jun 2015.
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
Scientific Reports, 50 (1):0 9981, 2015.
ISSN 2045-2322.
URL <https://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:0 022317, Aug 2016.
URL <https://link.aps.org/doi/10.1103/PhysRevA.94.022317>.

[Childs and Li, 2017]
Andrew M. Childs and Tongyang Li.
Efficient simulation of sparse markovian quantum dynamics.
Quantum Info. Comput., 170 (11–12):0
901–947, September 2017.
ISSN 1533-7146.
URL <https://dl.acm.org/doi/10.5555/3179568.3179569>.

[Cleve and Wang, 2016]
Richard Cleve and Chunhao Wang.
Efficient Quantum Algorithms for Simulating Lindblad Evolution.
arXiv e-prints, art. arXiv:1612.09512, December 2016.
URL <https://arxiv.org/abs/1612.09512>.

[Liu et al., 2020]
Zidu Liu, L. M. Duan, and Dong-Ling Deng.
Solving Quantum Master Equations with Deep Quantum Neural Networks.
arXiv e-prints, art. arXiv:2008.05488, August 2020.
URL <https://arxiv.org/abs/2008.05488>.

[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:0 022301, Jan 2000.
URL <https://link.aps.org/doi/10.1103/PhysRevA.61.022301>.

[Yoshioka et al., 2020]
Nobuyuki Yoshioka, Yuya O. Nakagawa, Kosuke Mitarai, and Keisuke Fujii.
Variational quantum algorithm for nonequilibrium steady states.
Phys. Rev. Research, 2:0 043289, Nov 2020.
URL <https://link.aps.org/doi/10.1103/PhysRevResearch.2.043289>.

[Harrow et al., 2009]
Aram W. Harrow, Avinatan Hassidim, and Seth Lloyd.
Quantum algorithm for linear systems of equations.
Phys. Rev. Lett., 103:0 150502, Oct 2009.
URL <https://link.aps.org/doi/10.1103/PhysRevLett.103.150502>.

[Knill et al., 2007]
Emanuel Knill, Gerardo Ortiz, and Rolando D. Somma.
Optimal quantum measurements of expectation values of observables.
Phys. Rev. A, 75:0 012328, Jan 2007.
URL <https://link.aps.org/doi/10.1103/PhysRevA.75.012328>.

[Lloyd, 1996]
Seth Lloyd.
Universal quantum simulators.
Science, 2730 (5278):0 1073–1078, 1996.
ISSN 0036-8075.
URL <https://science.sciencemag.org/content/273/5278/1073>.

[Nigro, 2019]
Davide Nigro.
On the uniqueness of the steady-state solution of the
Journal of Statistical Mechanics: Theory and Experiment,
20190 (4):0 043202, apr 2019.
URL <https://doi.org/10.1088

[Minganti et al., 2018]
Fabrizio Minganti, Alberto Biella, Nicola Bartolo, and Cristiano Ciuti.
Spectral theory of liouvillians for dissipative phase transitions.
Phys. Rev. A, 98:0 042118, Oct 2018.
URL <https://link.aps.org/doi/10.1103/PhysRevA.98.042118>.

[Nielsen and Chuang, 2016]
Michael A. Nielsen and Isaac L. Chuang.
Quantum Computation and Quantum Information (10th Anniversary
Cambridge University Press, 2016.
ISBN 978-1-10-700217-3.

[Barenco et al., 1995]
Adriano Barenco, Charles H. Bennett, Richard Cleve, David P. DiVincenzo, Norman
Margolus, Peter Shor, Tycho Sleator, John A. Smolin, and Harald Weinfurter.
Elementary gates for quantum computation.
Phys. Rev. A, 52:0 3457–3467, Nov 1995.
URL <https://link.aps.org/doi/10.1103/PhysRevA.52.3457>.

[Lacroix et al., 2020]
Nathan Lacroix, Christoph Hellings, Christian Kraglund Andersen, Agustin
Di Paolo, Ants Remm, Stefania Lazar, Sebastian Krinner, Graham J. Norris,
Mihai Gabureac, Johannes Heinsoo, Alexandre Blais, Christopher Eichler, and
Andreas Wallraff.
Improving the performance of deep quantum optimization algorithms
with continuous gate sets.
PRX Quantum, 1:0 110304, Oct 2020.
URL <https://link.aps.org/doi/10.1103/PRXQuantum.1.020304>.

[Nam et al., 2020]
Yunseong Nam, Yuan Su, and Dmitri Maslov.
Approximate quantum fourier transform with o(n log(n)) t gates.
npj Quantum Information, 60 (1):0 26, Mar
ISSN 2056-6387.
URL <https://doi.org/10.1038/s41534-020-0257-5>.

[Cleve> and Watrous, 2000]
R. Cleve and J. Watrous.
Fast parallel circuits for the quantum fourier transform.
In Proceedings 41st Annual Symposium on Foundations of Computer
Science, pages 526–536, 2000.

[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:0 011020, Feb 2021.
URL <https://link.aps.org/doi/10.1103/PhysRevX.11.011020>.

[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, 30 (3):0
1900052, 2020.

[Žnidarič, 2015]
Marko Žnidarič.
Relaxation times of dissipative many-body quantum systems.
Phys. Rev. E, 92:0 042143, Oct 2015.
URL <https://link.aps.org/doi/10.1103/PhysRevE.92.042143>.

[McClean et al., 2016]
Jarrod R McClean, Jonathan Romero, Ryan Babbush, and Alán Aspuru-Guzik.
The theory of variational hybrid quantum-classical algorithms.
New Journal of Physics, 180 (2):0 023023, feb
URL <https://doi.org/10.1088

[Peruzzo et al., 2014]
Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou,
Peter J. Love, Alán Aspuru-Guzik, and Jeremy L. O’Brien.
A variational eigenvalue solver on a photonic quantum processor.
Nature Communications, 50 (1):0 4213, 2014.
ISSN 2041-1723.
URL <https://doi.org/10.1038/ncomms5213>.

[Roggero and Baroni, 2020]
A. Roggero and A. Baroni.
Short-depth circuits for efficient expectation-value estimation.
Phys. Rev. A, 101:0 022328, Feb 2020.
URL <https://link.aps.org/doi/10.1103/PhysRevA.101.022328>.

[Rota et al., 2017]
R. Rota, F. Storme, N. Bartolo, R. Fazio, and C. Ciuti.
Critical behavior of dissipative two-dimensional spin lattices.
Phys. Rev. B, 95:0 134431, Apr 2017.
URL <https://link.aps.org/doi/10.1103/PhysRevB.95.134431>.

[Luo> et al., 2020]
Di Luo, Zhuo Chen, Juan Carrasquilla, and Bryan K. Clark.
Autoregressive Neural Network for Simulating Open Quantum Systems
via a Probabilistic Formulation.
arXiv e-prints, art. arXiv:2009.05580, September 2020.
URL <https://arxiv.org/abs/2009.05580>.

§ IMPLEMENTATION OF THE TRANSVERSE-FIELD ISING MODEL

We give here a specific circuit implementation of the controlled-$U^{2^j}$ operations in the case of the dissipative transverse-field Ising model, which is currently used as a benchmark for new numerical approaches to open quantum systems [Vicentini et al., 2019, Luo et al., 2020]. We assume that the computational basis encodes spin-1/2 eigenstates of the $Z$ operator. For $N$ spin, the model Hamiltonian is

\begin{equation}
\hat H = \frac{J}{4}\sum_{\langle j,k\rangle=1}^NZ_jZ_k+\frac{h}{2}\sum_{j=1}^NX_j\,,
\label{ising}
\end{equation}

where $\langle j,k\rangle$ denotes pairs of indices corresponding to first neighbor sites on the lattice being modeled. Dissipation is described by the jump operator $\hat A_j=\sigma^{(-)}_j=(X_j-iY_j)/2$. Both $\hat H$ and $\hat A_j$ are real-valued in the computational basis, so that $\hat H^T=H$, $\hat A_j^*=\hat A_j$, and $\hat A_j^\dagger=\hat A_j^T=\sigma^{(+)}_j=(X_j+iY_j)/2$.

For clarity, we omit the tensor products in Eq. (<ref>) from the notation in the formalism that follows. For this, we index the local operators on the left side of the tensor product with indices $j=1,\ldots,N$, and those on the right side with indices $j=N+1,\ldots,2N$. We attribute the index $j=0$ to the additional spin operators appearing on the left side of the tensor product in the definition of $M$, Eq. (<ref>). In this way, after simple algebra, the operator $M$ is rewritten as a Hamiltonian for an effective spin model with $2N+1$ spins

\begin{align}
M &= \frac{J}{4}\sum_{\langle j,k\rangle=1}^NY_0(Z_{j+N}Z_{k+N}-Z_jZ_k)+\frac{h}{2}\sum_{j=1}^NY_0(X_{j+N}-X_j)\nonumber\\
\label{Mising}
\end{align}

This effective Hamiltonian is $3$-local, i.e. it is written as $M=\sum_{\gamma}\alpha_{\gamma}M_{\gamma}$ where each $M_{\gamma}$ is a product of up to three Pauli operators.

As pointed out above, the controlled unitaries in the QPE part of the algorithm are implemented using the Suzuki-Trotter expansion. As a result, the implementation of the controlled unitary is reduced to the implementation of $C-e^{i\delta M_{\gamma}}$ gates, where $\delta$ is a small real-valued parameter. At the circuit level, this is achieved in terms of CNOTs, $\pm\pi/2$ single-qubit rotations along the $x$ or $y$ axes, and a controlled rotation by an arbitrary angle along the $z$ axis [Tacchino et al., 2020]. In Fig. <ref> we show the circuits for some selected two- and three-spin controlled operators, while those for the remaining terms occurring in the Suzuki-Trotter expansion can be similarly laid down.

@C=1em @R=1.0em @!R
(a)                 2

R_y(π/2)     1         1     R_y(-π/2)

@C=1em @R=1.0em @!R
(b)                 2

R_x(π/2)     1         1     R_x(-π/2)

R_y(π/2)         R_z(-2δ)         R_y(-π/2)

@C=1em @R=1.0em @!R
(c)                     3

R_x(π/2)     1                 1                 R_x(-π/2)

R_y(π/2)         1         1         R_y(-π/2)

R_x(π/2)             R_z(-2δ)             R_x(-π/2)

Quantum circuits implementing respectively (a) $C-e^{i\delta X\otimes Z}$, (b) $C-e^{i\delta Y\otimes X}$, and (c) $C-e^{i\delta Y\otimes X\otimes Y}$

Assuming to use a first-order Trotter expansion over $r$ Trotter steps, we have

\begin{equation}
U^j=e^{i2\pi j t_0 M}=\left(e^{i\frac{2\pi j t_0}{r}M}\right)^r\,,
\end{equation}

\begin{equation}
e^{i\frac{2\pi j t_0}{r}M}\simeq \prod_\gamma e^{i\frac{2\pi j t_0}{r}\alpha_\gamma M_\gamma}\,.
\label{Trotter1}
\end{equation}

This provides the final prescription to apply the $C-U^j$ as a sequence of gates of the form $C-e^{i\delta M_{\gamma}}$. Based on the circuits shown in Fig. <ref>, the gate cost for the implementation of the single Trotter step Eq. (<ref>) is $40N$ single-qubit gates, $42N$ CNOTs and 1 controlled $R_z$ gate. Twice as many gates would be needed if the second order Suzuki-Trotter expansion was employed instead of Eq. (<ref>).