A Perturbative Method for Nonequilibrium Steady State of Open Quantum Systems

# A Perturbative Method for Nonequilibrium Steady State of Open Quantum Systems

Tatsuro Yuge
Department of Physics, Osaka University, 1-1, Toyonaka, Osaka, 560-0043, Japan

Ayumu Sugita
Department of Applied Physics, Osaka City University, Sugimoto, Osaka, 558-8585, Japan
###### Abstract

We develop a method of calculating the nonequilibrium steady state (NESS) of an open quantum system that is weakly coupled to reservoirs in different equilibrium states. We describe the system using a Redfield-type quantum master equation (QME). We decompose the Redfield QME into a Lindblad-type QME and the remaining part . Regarding the steady state of the Lindblad QME as the unperturbed solution, we perform a perturbative calculation with respect to to obtain the NESS of the Redfield QME. The NESS thus determined is exact up to the first order in the system-reservoir coupling strength (pump/loss rate), which is the same as the order of validity of the QME. An advantage of the proposed method in numerical computation is its applicability to systems larger than those in methods of directly solving the original Redfield QME. We apply the method to a noninteracting fermion system to obtain an analytical expression of the NESS density matrix. We also numerically demonstrate the method in a nonequilibrium quantum spin chain.

## 1 Introduction

Establishing the statistical mechanics of the nonequilibrium steady state (NESS) is one of the most challenging problems in physics. In analogy to equilibrium statistical mechanics, a possible and desirable answer to this problem would be an ensemble theory with which we could write down the density matrix (or distribution function) for a NESS. To proceed in this direction, it is important to extract the characteristics of density matrices for NESSs by constructing and analyzing NESSs in various systems.

To realize a NESS in a concrete system, we should put the system in contact with heat baths or reservoirs that absorb the energy dissipated from the system. We thus have to treat open systems to analyze NESS. In quantum systems, one of the theoretical frameworks widely used for open systems is the quantum master equation (QME) [1], an equation of motion for the density matrix of the system. In fact, the QME is used in various fields of physics: e.g., quantum optics [1, 2], nuclear magnetic resonance [3], electron transfer in chemical physics and biophysics [4, 5], heat transport [6, 7, 8, 9, 10, 11, 12, 13, 14], electronic transport in mesoscopic conductors [15, 16, 17, 18, 19], spin transport [20], and nonequilibrium thermodynamics and statistical physics [21, 22, 23, 24, 25]. Therefore, the QME is a reliable approach to investigating NESS in various systems (c.f. Ref. [26]). In some solvable models, analytical expressions of the NESS of the QME are obtained [27, 28, 29, 30, 31, 32].

Apart from solvable models, however, there appears a difficulty in analyzing the NESS of the QME. This is because one should treat a huge number of components in the QME; the number of elements of a density matrix (mixed state) is the square of that of a wave function (pure state). This restricts available system size to be relatively small when one employs direct methods (direct diagonalization or time integration) in solving the QME, even though efficient numerical methods are developed [4, 33].

A method of treating larger systems in the framework of the QME is the stochastic unraveling of QME (also known as quantum jump, quantum trajectory, and Monte Carlo wave function method) [1, 2, 34, 35, 36, 37, 38, 39, 41, 42, 40, 5, 10, 11]. In this method, since one treats a wave function instead of a density matrix, one can investigate systems larger than those in the direct methods. This method, however, requires a Monte Carlo calculation; i.e., one must take an average of many times of iterative computations. There are also other methods. In Ref. [13], a density matrix renormalization group method with a matrix product operator ansatz is proposed. This method assumes the (rather phenomenological) local action of the dissipators associated with the baths. In Ref. [43], a many-particle Green function method is proposed. This requires the approximation of weak internal coupling or small system size.

In this paper, we propose another method of calculating the NESS of the QME. We use a perturbation theory with respect to the coupling strength (pump/loss rate) between the system and the reservoirs. Our method assumes only the weak system-reservoir coupling, which is the same assumption as that of the QME. In this method, the unperturbed part of the NESS has only diagonal elements in the energy eigenstate basis. Thus, the number of elements of the unperturbed part is on the same order as that of the wave function. We can therefore treat systems as large as those in a stochastic unraveling method. Moreover, since this method does not require Monte Carlo calculation, it is expected to be faster than the unraveling method. We also note that the Markovian approximation use in the QME holds exactly for the steady-state solution, since the approximation is good when the system changes slowly [22].

In the present paper, we describe the method in a setup of a system that is weakly coupled to two reservoirs. In the next section, we explain the setup and QME. In Sect. 3, we give the main result, i.e., the perturbative solution for the NESS of the QME. We also show an advantage of our method in numerical computation. Furthermore, we explain how we calculate currents in this method. In Sect. 4, we have two examples. One is a noninteracting fermion system, where we derive an analytical expression of the NESS by the method. The other is a numerical computation of the NESS in a quantum spin chain, where we demonstrate the validity of the method. We devote Sect. 5 to the concluding remarks.

## 2 Setup

We consider transport phenomena (energy transport, particle transport, etc.) in the quantum system S. S is in contact with two reservoirs (heat baths), L and R (see Fig. 1), and is thus an open system. We assume that the dimension of the Hilbert space associated with S is finite (). Each reservoir is sufficiently large compared with S and is in an equilibrium state characterized by its own inverse temperature , chemical potential , and so on (). We also assume that the coupling between the system S and each reservoir is weak. The Hamiltonian of the total system (L+S+R) reads

 ^Htot=^HS+∑b=L,R(^Hb+u^HSb). (1)

The first term is the Hamiltonian of S, which may be degenerate. We denote the eigenenergy as and the corresponding eigenstate as , where is the label for distinguishing the degenerate states. The second term is the Hamiltonian of the reservoir . The third term is the system-reservoir coupling Hamiltonian. The small prefactor represents the weak system-reservoir coupling assumption.

### 2.1 Quantum master equation

Here, we explain the quantum master equation (QME) to fix the notation in the present paper. The starting point is the Liouville-von Neumann equation of the total system:

 d^ρtot(t)dt=1iℏ[^Htot,^ρtot(t)], (2)

where represents the density matrix of the total system. We can derive the QME, the equation of motion of the system S, by applying the Born-Markov approximation [1] to Eq. (2). The QME in the Schrödinger picture reads

 d^ρ(t)dt=L^ρ(t). (3)

Here, is the reduced density matrix of S, and is the trace over the reservoir . The superoperator (QME generator) is given by

 L≡L0+v∑bLb, (4)

where is a parameter that controls the pump/loss rate due to the coupling to the reservoirs, describes the unitary time evolution, and

 Lb^ρ≡−1ℏ2∫∞0dt′Trb[^HSb,[˘HSb(−t′),^ρ⊗^ρb]] (5)

is the dissipative part associated with the reservoir . Here, with is an operator in the interaction picture, and is the equilibrium state of the reservoir . The explicit form of depends on situations; the canonical ensemble or the grand canonical ensemble is usually used. In deriving these equations, we assumed . Equation (3) with Eq. (4) is called Redfield QME. We note that the QME is valid up to first order in () since the Born approximation is a second-order approximation with respect to .

We now rewrite the QME into a more tractable form. To this end, we first assume that the system-reservoir coupling Hamiltonian has the form

 ^HSb=∑λ^Xb,λ⊗^Yb,λ, (6)

where is a self-adjoint operator of S that is defined locally in the regions (hatched in Fig. 1) near the reservoir , and is a self-adjoint operator of the reservoir . We next decompose into eigenoperators of [1]:

 ^Xb,λ =∑ω^X(ω)b,λ, (7) ^X(ω)b,λ ≡∑i^Π(Ei)^Xb,λ^Π(Ei+ℏω), (8)

where runs over the eigenenergy differences, and is the projection operator onto the eigenspace with the eigenenergy . Substituting Eqs. (6) and (7) into Eq. (5), we have

 Lb^ρ =1ℏ2∑λ,ν∑ωΞbλν(ω){^X(ω)†b,ν^ρ^Xb,λ−^Xb,λ^X(ω)†b,ν^ρ}+h.c. =1ℏ2∑λ,ν∑ω,ω′Ξbλν(ω){^X(ω)†b,ν^ρ^X(ω′)b,λ−^X(ω′)b,λ^X(ω)†b,ν^ρ}+h.c., (9)

where is the Fourier-Laplace transform of the reservoir correlation function:

 Ξbλν(ω)≡∫∞0dte−iωtTrb[˘Yb,λ(t)^Yb,ν^ρb]. (10)

In Eq. (9), we used and .

### 2.2 Decomposition of QME

From Eq. (9), we can decompose into the part with and the remaining one: , where

 LSAb^ρ ≡1ℏ2∑λ,ν∑ωΞbλν(ω){^X(ω)†b,ν^ρ^X(ω)b,λ−^X(ω)b,λ^X(ω)†b,ν^ρ}+h.c. (11)

Here, is a Lamb shift Hamiltonian, , , and the overlines stand for the operation of taking the complex conjugate. Since is positive definite, gives a Lindblad-type QME. That is, the -omitted superoperator is a generator of a completely positive dynamical semigroup [44, 45]. The approximation of omitting the terms is called secular approximation (SA) or rotating wave approximation [1]. The SA is an approximation that extracts a slow part of the dynamics; one can obtain by averaging in time in the interaction picture or by omitting the fast oscillating terms in in the interaction picture. However, it is known that the internal current vanishes in the NESS of the SA-QME [9]. Therefore, the SA-QME is not appropriate for analyzing the steady state itself (however, note Ref. [46]), and one often uses any of the original Redfield QME, alternatively approximated Lindblad QME [9], or axiomatic Lindblad QME [44, 45] for nonequilibrium situations.

In the present study, we use the Redfield QME. From the above argument, we have a decomposition of the QME generator:

 L=LSA+vR, (12)

with . In the next section, we develop a perturbative method of calculating the NESS of the Redfield QME (3), where we regard and respectively as the unperturbed and perturbation parts.

### 2.3 Liouville space and projection superoperator

To summarize the points so far, we have introduced the Redfield QME (3) that describes the dynamics of the open quantum system S. The QME generator is written as Eq. (4) or (12). Before going to the main result in the next section, here, we make two more preliminaries.

One is the Liouville space , which is the set of all the linear operators on . Since the dimension of is finite, any in is a trace class operator. We can define the Hilbert-Schmidt inner product in as for any , where is the trace in . With this inner product, is a Hilbert space. We refer to any linear operator on as a superoperator (which includes , and so on). We note that the dimension of is given as . We define the adjoint of a superoperator such that holds for any .

The other is the following projection superoperator , which is defined by

 (13)

In the matrix representation of any in the basis of the energy eigenstates of , extracts only the matrix elements constructed from the eigenstates with the same eigenenergies. We also define the projection superoperator . By using these superoperators, we decompose into the subspace and its orthogonal complement . We note that is on the same order as (), particularly, if eigenenergies of are not degenerate.

By using the above projection superoperators, we can show the following relations:

 PL0 =L0P=0, (14) PLSAbQ =QLSAbP=PLSAQ=QLSAP=0, (15) PRbP =PRP=0. (16)

Equation (15) leads to , which implies that the eigenvectors of are classified into those belonging to and . Equation (16) implies that holds, and that maps operators in to ones in ().

## 3 Perturbative Method for NESS

In the present paper, we are interested in the NESS of the QME. By substituting into the QME (3), we obtain the steady-state equation (since is time-independent),

 L^ρss=0. (17)

Thus, our task is to solve the zero-eigenvalue problem of .

For this purpose, we take the following strategy that consists of two steps. In the first step, we solve the zero-eigenvalue problem of , i.e., . Here, we assume that there exists a unique . See Refs. [21, 47] and [48] for the conditions for this assumption. We may solve this problem by either numerical or analytical methods. As mentioned earlier, the steady state in the secular approximation is not suitable for analyzing NESS. However, combined with the following second step, is a useful starting point to obtain NESS.

In the second step, taking as the unperturbed eigenvector, we perform the first-order perturbation calculation with respect to [see Eq. (12)]. Then, we have an correction to , which gives a perturbative solution of Eq. (17). This perturbation theory is valid up to the first order in , which is the same order as that in the QME (3).

This strategy shares the spirit in Ref. [26] of solving the steady-state equation [Eq. (17)] by a perturbative method that is valid up to . The difference lies in the decomposition into unperturbed and perturbation parts.

### 3.1 Perturbation theory with respect to vR

Now, we analyze in detail the perturbation theory mentioned in the second step above. First, we define the eigenvalues and the corresponding left and right eigenvectors and of the unperturbed generator as

 L†SA^ℓk =¯¯¯¯¯λk^ℓk, (18) LSA^rk =λk^rk. (19)

We assign to the zero eigenvalue; i.e., , (identity operator on ), and is the steady solution of .

Next, we apply a formal perturbation theory to the zero eigenvalue and eigenvectors of with respect to . Note that is nondegenerate because we assume that is uniquely determined. Therefore, similarly to the case of quantum mechanics [49], we can use the perturbation theory for the nondegenerate case to obtain the first-order terms , , and (corrections to , , and , respectively) as

 δλ =vTrS[^ℓ†0R^r0], (20) ˆδℓ =v∑k≠0¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯⎛⎜⎝TrS[^ℓ†0R^rk]λ0−λk⎞⎟⎠^ℓk, (21) ˆδr =v∑k≠0TrS[^ℓ†kR^r0]λ0−λk^rk. (22)

Since , Eq. (20) and the numerators in Eq. (21) have the form of . By using and the trace-preserving property of the Redfield and Lindblad QMEs, we have . Therefore, the correction terms and vanish. This fact indicates that the eigenvalue and left eigenvector remain respectively zero and in this perturbation theory, and thus the corresponding right eigenvector is indeed the steady state of .

Next, we investigate the correction term to the right eigenvector. As mentioned earlier concerning Eq. (15), each eigenvector of belongs to either or . In particular, belongs to since ( is a density matrix). This fact with Eq. (16) leads to (). Therefore, the terms with in Eq. (22) vanish. Thus, we obtain

 ˆδr=−v∑k∈QTrS[^ℓ†kR^ρSAss]λk^rk, (23)

where the sum runs over the labels whose eigenvectors belong to . To rewrite further, we evaluate , , and in the above equation. To this end, we note that it is sufficient to evaluate them to because includes in the prefactor. Therefore, we can replace , , and with the eigenvalue and eigenvectors of . Since the eigenvalue is and both the left and right eigenvectors are , we finally obtain

 ˆδr =−iℏv∑i≠j∑mi,mj⟨Ei,mi|R^ρSAss|Ej,mj⟩Ei−Ej|Ei,mi⟩⟨Ej,mj|. (24)

This result implies two important points. One is that all the denominators in Eq. (24) are . This is necessary for the validity of this perturbation theory since if one of the denominators is on the order of , becomes on the order of , which is contradictory to the assumption that this is the first-order term. The other point is that, to calculate , we have to determine not all the eigenvalues and eigenvectors of but only the right eigenvector for the zero eigenvalue.

We thus construct a perturbative method of calculating the NESS of the QME (3):

 ^ρss=^ρSAss+ˆδr (25)

(or equivalently, and ) with Eq. (24). This is the main result in the present paper. We note that this result is very similar to the result in Ref. [50].

### 3.2 Advantage in numerical computation

The method of Eq. (25) requires information on , which is the right eigenvector of () for the zero eigenvalue. Therefore, when we use this method in calculating the NESS of the QME (3), we need not directly solve the full steady-state equation (zero-eigenvalue problem) of , but we have to solve the reduced zero-eigenvalue problem of . Since the dimension of the matrix () is much smaller than that of (), this method has an advantage in the numerical computation of .

The computation of with Eq. (24) requires information on . Although the matrix () is not significantly smaller than , the computational cost of the matrix-vector multiplication is much smaller than that of solving the eigenvalue problem.

In Ref. [9], Wichterich et al. proposed a Lindblad-type QME that has nonvanishing internal current in NESS, where they assumed weak internal couplings and high temperature. In contrast, our method does not require these assumptions. When using a stochastic unraveling method [1, 2, 34], one can reduce the problem size from to . However, it requires a Monte Carlo calculation; i.e., one must take the average of many runs of computation. By comparison, our method determines the NESS in a single run.

### 3.3 Current

Here, we consider energy currents in the present formulation. As shown in Fig. 2, we virtually divide the system S into two parts; the left part includes the interaction region with the reservoir L, while the right part includes the interaction region with the reservoir R. Correspondingly to this division, we can write the system Hamiltonian as , where and are respectively the Hamiltonian of the left and right parts. The interaction between the left and right parts can be included in either or . We then define the energy current from the left part to the right part as the unit-time loss of the left energy due to the interaction with the right part: . This is equivalent to the unit-time gain of the right energy due to the interaction with the left part, , because the energy conservation holds if the interactions with the reservoirs are absent ().

We note that holds (because and ), which indicates that the internal current vanishes in the SA-QME, as mentioned earlier. Using this fact, we have the steady-state average of the current in the Redfield QME as

 ⟨^Jl→rE⟩ss =v∑i≠j∑mi,mj⟨Ej,mj|^Hl|Ei,mi⟩⟨Ei,mi|R^ρSAss|Ej,mj⟩ =vTr[^HlR^ρSAss]. (26)

Here, we used Eq. (24) for in the second line and in the third line.

We next investigate the energy current from the reservoir L to the system S. We define the current operator as the unit-time energy gain of S due to the coupling with the reservoir L: . We can divide it as , where is the current from the reservoir L to the left part of S, and is the current from the reservoir L to the right part (see Fig. 2). However, holds if the right part of S is sufficiently separate from the interaction region with the reservoir L (the hatched region in the left part in Fig. 2). We note that this condition is usually satisfied (except for extremely small systems). We can show the operator identity as

 ^JL→rE =v2ℏ2∑λ,ν∑ωΞLλν(ω){^XL,λ^Hr^X(ω)†L,ν−^Hr^XL,λ^X(ω)†L,ν}+h.c. =0. (27)

Here, we used the commutability of with (because of the lack of overlap between the right part and the interaction region with the reservoir L).

Similarly, we can define the energy current from the system S to the reservoir R , and show that it is equivalent to the current from the right part to the reservoir R because the current from the left part to the reservoir R is absent ().

By using these facts, we can show the continuity equations of the energy in the Heisenberg picture: and . In the steady state, we thus obtain as expected.

## 4 Examples

### 4.1 Analytical expression of NESS in noninteracting fermion system

We consider a noninteracting fermion system that is connected with two reservoirs, as an example. The reservoirs are in equilibrium states at different temperatures and chemical potentials, so that the particle and energy currents flow in the system. Using the method of Eqs. (25) and (24), we derive an analytical expression of the density matrix of the NESS.

We consider a system S of fermions that move on a one-dimensional -site lattice:

 ^HS=N∑l=1εl^d†l^dl+N−1∑l=1(thopl^d†l^dl+1+h.c.), (28)

where is the energy level of the th site, is the transfer probability amplitude between the th and th sites, and and are respectively the creation and annihilation operators of the fermion at the th site. We can diagonalize this Hamiltonian by an appropriate linear transformation:

 ^HS=N∑k=1ℏωk^c†k^ck, (29)

where

 ^dl=N∑k=1Wlk^ck,  ^d†l=N∑k=1¯¯¯¯¯¯Wlk^c†k, (30)

and is the transformation matrix. For simplicity, we assume that if . In this case, is nondegenerate.

At the left () and right () ends of the lattice, the system S is respectively coupled to the reservoirs L (at the inverse temperature and the chemical potential ) and R (at the inverse temperature and the chemical potential ). We assume that these reservoirs are free fermion systems:

 ^Hb=∑qεb,q^c†b,q^cb,q, (31)

where and are the creation and annihilation fermion operators in the reservoir , respectively. The coupling Hamiltonian is given by

 ^HSb =∑q(ℏξb,q^d†l(b)^cb,q+h.c.) =∑k,q(ℏζkbq^c†k^cb,q+h.c.). (32)

In the second line, we used Eq. (30). Here, is the coupling constant, , if , and if . We note that, unlike Eq. (6), is not self-adjoint. However, the generic formulation in Sects. 2 and 3 is unmodified by this difference. In this case, it is convenient to use a slightly different definition of eigenoperators:

 ^c(ω)k ≡∑i^Π(Ei)^ck^Π(Ei+ℏω), ^c†(ω)k ≡∑i^Π(Ei)^c†k^Π(Ei−ℏω). (33)

In this definition, holds. In particular, in the case of the noninteracting fermion system, and . In almost the same manner as in Sect. 2, we obtain the Redfield QME generator , where with and

 Lb^ρ=12∑k1,k2∑ωΓb(ω) {¯¯¯¯¯¯Wl(b)k1Wl(b)k2f+b(ω)(^c†(ω)k1^ρ^ck2+^c†k1^ρ^c(ω)k2−^ck2^c†(ω)k1^ρ−^ρ^c(ω)k2^c†k1) +Wl(b)k1¯¯¯¯¯¯Wl(b)k2f−b(ω)(^c(ω)k1^ρ^c†k2+^ck1^ρ^c†(ω)k2−^c†k2^c(ω)k1^ρ−^ρ^c†(ω)k2^ck1)}.

Here, is the reservoir spectral function and . In the above equation, we omitted terms proportional to the imaginary part of [7], since they do not have a significant contribution. Also, the SA-QME generator is given as , where

 LSAb^ρ=12∑kγkb{f+b(ωk)(2^c†k^ρ^ck−^ck^c†k^ρ−^ρ^ck^c†k)+f−b(ωk)(2^ck^ρ^c†k−^c†k^ck^ρ−^ρ^c†k^ck)},

with .

We can solve the zero-eigenvalue problem of the above with the corresponding eigenvector of the form of , where

 ^ρk =1∑bγkb∑bγkb[f−b(ωk)^ck^c†k+f+b(ωk)^c†k^ck]. (34)

This is equivalent to a weighted average of the equilibrium states with the parameters corresponding to the reservoir L or R: . (Note that the equilibrium state for this model at the inverse temperature and the chemical potential is written as ). We can also express as

 ^ρSAss =1Zssexp[−∑kGk^c†k^ck], (35)

where . This result is equivalent to the result in Sect. III A of Ref. [51].

By using Eq. (24) with this result, we obtain the first-order correction of the NESS density matrix:

 ˆδr= vi2∑k1≠k21ωk1−ωk2[{F+k1k2^ck1^c†k2^ρSAss−F−k1k2^ρSAss^ck1^c†k2}−h.c.], (36)

where

 F±k1k2≡∑bγk2b∑bγk2bf±b(ωk2)∑b¯¯¯¯¯¯Wl(b)k1Wl(b)k2Γb(ωk1)f±b(ωk1).

Equations (35) and (36) are respectively analytical expressions of the diagonal and off-diagonal elements of the NESS density matrix in the noninteracting fermion system.

### 4.2 Numerical calculation in quantum spin chain

As another example, we apply our method to the numerical computation of the NESS in a one-dimensional Ising spin chain subject to a magnetic field. The Hamiltonian of the chain system S reads

 ^HS=−JN−1∑l=1^σzl^σzl+1−hxN∑l=1^σxl−hzN∑l=1^σzl, (37)

where () is the -component of the Pauli matrices at the th site, is the Ising coupling constant, is the -component of the magnetic field, and is the total site number in the chain. At the left () and right () ends of the chain, the system S is coupled to the heat reservoirs L (at the inverse temperature ) and R (at the inverse temperature ), respectively. We assume that these reservoirs are composed of free boson particles:

 ^Hb=∑kεb,k^a