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

Simulating the dynamics and the non-equilibrium steady state of an open quantum systems 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 relaxation times or strong quantum correlations in the transient states. Here, we develop an efficient quantum algorithm for the direct calculation of the non-equilibrium steady state on a fault-tolerant quantum computer. The algorithm makes use of quantum phase estimation to evaluate 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, using existing methods for the quantum measurement of observables. 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 [1][2][3]. 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 [4,5] to the accurate modeling of noise in modern quantum computing platforms [6]. The numerical simulation of the dynamics of a manybody 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 [7] to the solution of the Lindblad-Von-Neumann master equation for the density matrix [8] -describing the case of a Markovian environment -include in addition to the exact numerical solution of the resulting differential equations, tensornetwork methods [9][10][11], projector [12] and time-dependent variational [13][14][15][16] quantum Monte Carlo, various quantum trajectory approaches, as well as methods introducing various levels of approximation [17][18][19][20][21]. 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 [22][23][24], and specific real-space decimation schemes [25]. 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 integrating the time-dependent master equation for the density matrix have recently been proposed [26][27][28][29][30][31][32][33][34][35][36]. 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 scheme. While some special cases of system-environment interaction have been shown to translate into a unitary stochastic quantum evolution [26], more generally the dynamics could be addressed either by an effective dynamics over an augmented Hilbert space, or by suitably approximating the non-unitary evolution with a unitary transformation over short times. Similarly, a time evolution quantum algorithm was proposed to compute the thermal equilibrium state of a system in interaction with a thermal bath [37]. Recently, a hybrid quantum algorithm has been proposed for the direct calculation of the NESS, based on a variational ansatz [38].
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. Here, we propose a different approach consisting in the direct quantum computation of the NESS without integrating the system dynamics, by directly solving the linear system of equations obeyed by the steady state. The algorithm uses quantum phase estimation (QPE) -similarly to the HHL algorithm [39] for the solution of linear systems of equations -and leverages the known spectral properties of the generator of the dissipative dynamics of the system. It then provides an approximation of the NESS density matrix. It is shown that in cases in which the Liouvillian gap -i.e. the asymptotic relaxation rate towards the NESSdecreases 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 [40], is it finally shown how to directly measure 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 and runtime is performed. In Section 4 the formalism to estimate expectation values on the computed NESS is derived. Section 5 contains a numerical analysis of the algorithm applied to a simple driven-dissipative problem. The conclusions are presented in Section 6.

Formalism
The dynamics of an open quantum system interacting with a Markovian environment is described by the Lindblad-Von Neumann master equationρ = L(ρ) [8], whereρ is the density operator, and the liouvillian super-operator is expressed as (setting = 1) Here,Ĥ is the system Hamiltonian andÂ 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Ĥ and theÂ 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 [41]. The case of global jump operators, describing transitions between eigenstates ofĤ, 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 L(ρ 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 [42], as it will always be the case on a quantum circuit implementation of the problem. We adopt the Choi isomorphism (also known as vectorization) to map operators on the Hilbert space H onto vectors on the tensor product space H ⊗ H. Given an orthonormal basis, a density operatorρ = jk ρ jk |j k| maps onto the vector |ρ = jk ρ jk |j ⊗ |k (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 H ⊗ H, defined as Here, we will assume that the Hilbert space H of the problem has dimension 2 N . Then the vector |ρ 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 λ j are in general complex, and the nonzero eigenvalues are characterized by Re(λ j ) < 0. The right and left eigenvectors associated to the zero-eigenvalue (here assumed unique) are respectively the NESSρ ss and the identity matrix, i.e. L|ρ ss = 0 and 1|L = 0 [8,43]. The latter relation implies also L † |1 = 0. We define a Hermitian operator in a 2 2N +1 -dimensional space as From the spectral properties of the Liouvillian, we gather that the operator M has two eigenvectors with zero eigenvalue, namely |η 0 = |0 |1 and |η 1 = |1 |ρ ss . From the assumption that the NESS is unique, we can conclude that these will be the only eigenvectors with zero eigenvalue of M . More generally, we denote with |η j and φ j the eigenvectors and the corresponding eigenvalues of M . Given the structure of M , for each eigenvalue φ j associated to the eigenvector |η j = |0 |a j +|1 |b j , the vector |η j = |0 |a j −|1 |b j is also eigenvector of M with eigenvalue −φ j . The eigenvector condition also implies the relations L † L|b j = ϕ 2 j |b j and LL † |a j = ϕ 2 j |a j . Hence, the eigenvalues of M are related to the eigenvalues λ j of the Liouvillian L through the relation ϕ 2 j = |λ j | 2 . This link shows that the Liouvillian gap g, i.e. the minimal value of |Re(λ j )| among all nonzero eigenvalues, is also a measure of 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 where X and Y are Pauli operators. If L H and L A can be efficiently encoded in terms of quasi-local operators, then M can also be efficiently encoded. From Eq. (2), 3 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 of the NESS solver is illustrated in Fig. 1. It comprises an input-state preparation stage, a QPE stage, and a measurement stage.
· · · · · · · · ·  Figure 2: Quantum circuit P for the preparation of the initial state |ξ for the second register. The circuit takes as input the (2N + 1)-qubit state |0 and outputs the state |ξ in Eq. (6). Controlled-Hadamard gates are implemented following the general procedure for arbitrary controlled gates [44,45], and controlled-arbitrary gates are available within current superconducting circuit technology [46].
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 [44], QPE takes as an input the state |0 |η j , where |η j is an eigenstate of a unitary operator U with eigenvalue e 2πiϕ j , assuming ϕ j ∈ [0, 1]. If 2 t ϕ j is a t-bit integer, then QPE outputs the state |ϕ j |η j . If 2 t ϕ j is not an integer, then QPE will produce with high probability the state |φ j |η j , where 2 tφ j is the best t-bit integer estimate of 2 t ϕ j from below.
For the present case, we set U = e 2πiM t 0 , where the real parameter t 0 is set so that the whole spectrum of M t 0 is included in the [0, 1] interval. The algorithm assumes that black-box circuit elements are available to efficiently perform controlled-U 2 j operations. Given the quasi-local character of M in the computational basis, U can be encoded using algorithms to simulate the time evolution of sparse Hamiltonians [41,47,48], and controlled-U 2 j operations are implemented through standard techniques [44].
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 |0 |ξ of the QPE stage of the algorithm must be chosen in such a way that |ξ has a significant overlap with the zero-eigenspace of M spanned by |η 0 and |η 1 . Here we set 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 |1 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 (6) when input with the initial state |0 , is shown in Fig. 2. Controlled Hadamard gates can be implemented with a limited number of gates following the prescriptions for arbitrary controlled operations [44,45], and controlled-arbitrary gates are available within current superconducting circuit technology [46]. Expanding |ξ on the basis of the eigenvectors |η j of M , the state of the full circuit following the preparation stage is where we have indicated explicitly the expansion of the state |1 |0 on the basis of eigenstates |η j of M , and accounted for the fact that states |1 |0 and |η 0 = |0 |1 are orthogonal. The first term of the sum in the last line is the one that will be untouched by the QPE, as both |η 0 and |η 1 are eigenvectors of M associated to the zero eigenvalue. In general, the other eigenvalues of M will be such that 2 t ϕ j is not a t-bit integer. Then, the QPE stage of the circuit applied to the state |0 |η j will output the state k |k |η j , so that the state of the circuit at the output of the QPE stage is where After the QPE stage, the first register is measured in the computational basis, and the algorithm is successful if the eigenvalue 0 is measured. Then, the output state is (up to a normalization constant of order 1) where c 1 ∈ 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 (10) is the result being sought, while the remaining term is an unwanted error proportional to the eigenstates |η j of M other than those with zero eigenvalues. We will show in the next Section that the probability amplitude associated to this spurious term decays as O(2 −t ) and can thus be made arbitrarily small by increasing the number of qubits t.

Success probability and runtime
From (8), the probability of measuring zero on the first register is 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 (10), it is however useful to study the behavior of |α (j) 0 | 2 . Simple algebra leads to |α (j) The function (12) is plotted in Fig. 3(a). As already pointed out, the value of t 0 is set in order to have ϕ j < 1. The value of |α (j) 0 | 2 can be large only when ϕ j approaches an integer value, and the value of t 0 can be set such that the only relevant case is when ϕ j ∼ 0. The peak at φ j ∼ 0 has width O(2 −t ) and the envelope of the function in the vicinity of ϕ j ∼ 0 decays as 2 −2t , implying that the probability |α (j) 0 | 2 is significantly nonzero only for eigenvalues ϕ j 2 −t . Then, the number of quantum bits t must be chosen such that the Liouvillian gap g, defined here as the eigenvalue ϕ j with the smallest nonzero absolute value, fulfills 2 −t g. The algorithm can thus distinguish between the zero and the closest nonzero eigenvalues provided t log(1/g).
The inverse gap for most model Hamiltonians and couplings to the environment scales polynomially with the system size, except at critical points where it can scale exponentially [49]. In the case of polynomial scaling, choosing a number of qubits t ∼ O(log(N )) is therefore sufficient to estimate the NESS with high accuracy. To estimate how the runtime of the algorithm scales with N , we first observe that both the initial state preparation P and the Hadamard stage H ⊗t have a number of gates that scales as O(N ), while the QFT stage has O(N 2 ) gates. The remaining factor determining the runtime is the application of the controlled-U operations, which amounts to the application of a U 2 t = e 2πi2 t M t 0 to the (2N + 1)-qubit register. The application of U 2 t can be executed in time O(2 t s 2 ) where s is the sparsity index of M [47]. For systems with short-range two-body interactions s = O(1), while long-range interactions will result in s = O(N ). Hence, for an N -body open quantum system with two-body interactions and power-law scaling of the inverse Liouvillian gap, the overall computational weight of the present algorithm will scale as O(2 t s 2 ) ∼ O(N ) for short-range interactions and at most as O(N 3 ), 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 [44] or a variational quantum eigensolver [50,51].

Estimate of expectation values
Once the first register has been measured in the state |0 , it is possible to estimate quantum mechanical expectation values from the density matrixρ ss . For a given observableÔ, the expectation value is computed as Ô = Tr(Ôρ ss )/Tr(ρ ss ). As before, we assume that the Hermitian operatorÔ 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 |(Ôρ) = 1 ⊗Ô|ρ , and the super-operator O = 1 ⊗Ô is also Hermitian. The trace average is then expressed as a matrix element Tr(Ôρ ss ) = 1|O|ρ ss = ρ ss |O|1 . (13) To take advantage of the representation of the NESS at the output of the algorithm, we define the Hermitian operator acting on the (2N + 1)-qubit second register. Then, while, due to the structure of the operator Q, we have η 0 |Q|η 0 = η 1 |Q|η 1 = 0.
The trace average can then be estimated as the expectation value of 1 ⊗ Q on the output state |ψ 3 where we used the facts that the expectation value of O is real-valued, η 0 |Q|η 1 = η 1 |Q|η 0 , and c 1 ∈ R. In particular, by settingÔ = I, the measurement returns the quantity c 1 Tr(ρ ss ), which can be then used to eliminate the factor c 1 and obtain the expectation value Ô = Tr(Ôρ ss )/Tr(ρ ss ). The second and third term on the right hand side of Eq. (16) are the spurious terms resulting from the finite number of qubits t in the QPE. They are respectively of order O(|α . 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. (13) and (16) prove that the problem of computing the trace average ofÔ maps onto that of estimating the expectation values of O and 1 onto the output quantum state |η ss of the present algorithm. This estimation is achieved through methods for the quantum measurement of expectation values [40] and requires, similarly to QPE, the possibility to efficiently encode the unitary operator W = e 2πiQt 0 which again results from the assumptions on the observableÔ.

Results
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Ĥ = hσ x and the density matrix obeys the Lindblad-Von Neumann master equatioṅ 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. 1 and 2 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ρ Results are displayed in Fig. 3(b). The accuracy of the quantum algorithm scales as 2 −2t , as expected from Eqs. (10) and (12). In Fig. 4 the expectation values ofσ y andσ z , estimated according to Eq. (16), 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 α (j) 0 in Eq. (16) results in a relative error scaling as δ ∼ 2 −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 solution of the NESS problem for 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 expectation value of an arbitrary observable onto the NESS by using the established method for the quantum measurement of expectation values. 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. With the availability of scalable fault-tolerant quantum hardware, this 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.