Local master equations may fail to describe dissipative critical behavior

Local quantum master equations provide a simple description of interacting subsystems coupled to different reservoirs. They have been widely used to study nonequilibrium critical phenomena in open quantum systems. We here investigate the validity of such a local approach by analyzing a paradigmatic system made of two harmonic oscillators each in contact with a heat bath. We evaluate the steady-state mean occupation number for varying temperature differences and find that local master equations generally fail to reproduce the results of an exact quantum-Langevin-equation description. We relate this property to the inability of the local scheme to properly characterize intersystem correlations, which we quantify with the help of the quantum mutual information.


I. INTRODUCTION
Quantum master equations have been instrumental in the study of open quantum systems since their introduction by Wolfgang Pauli in 1928 [1]. They offer powerful, yet approximate, means to describe the time evolution of the reduced density operator of quantum systems coupled to external environments [2][3][4][5]. They allow the analysis of the dynamics of both diagonal density matrix elements (populations), involved in thermalization processes, and of nondiagonal density matrix elements (coherences), associated with dephasing phenomena. As a consequence, they have found widespread application in many different areas, ranging from quantum optics [6] and condensed matter physics [7] to nonequilibrium statistical mechanics [8] and quantum information theory [9].
However, the form of the quantum master equations employed in these studies is often postulated. Their validity is thus not completely clear a priori. This is especially true for boundary-driven processes where the system of interest is coupled to several reservoirs. In this case, it has recently been shown that local master equations, that are commonly used to examine nonequilibrium phase transitions [17][18][19][20][21][22][23][24][25][26][27][28], may violate the second law of thermodynamics [29] and give rise to nonphysical results, such as incorrect steady-state distributions or nonzero currents for vanishing bath interactions [30][31][32][33][34][35][36][37][38][39][40], even in the limit of small bath couplings. These inconsistencies are related to the fact that local quantum master equations, whose total dissipator is simply λ T 1 T 2 a 1 a 2 FIG. 1. Coupled-oscillator model. Two quantum harmonic oscillators interact with each other with interaction strength λ. Each of them is weakly coupled with a heat bath with respective temperature Tj, (j = 1, 2). A nonequilibrium steady state is established when the two temperatures are different and heat flows from one oscillator to the other.
the sum of the single-bath dissipators, incorrectly neglect bath-bath correlations, which are induced by intersystem interactions, in contrast to global quantum master equations [29][30][31][32][33][34][35][36][37][38][39][40]. Interestingly, the local approach has been shown to provide a better description of quantum heat engines than the global approach in some parameter regimes [35]. Meanwhile, the validity of Lindblad quantum master equations has, for example, been discussed in the context of quantum transport [41,42], quantum relaxation [43,44], and entanglement generation [45]. But these results cannot be straightforwardly extended to nonequilibrium phase transitions as the considered models do not exhibit critical behavior.
In this paper, we examine the accuracy of a quantummaster-equation description of dissipative critical phenomena by analyzing an exemplary system consisting of two interacting harmonic oscillators, each weakly coupled to a thermal reservoir. This system naturally appears in many areas, most notably in cavity optomechanics [46]. Many-body superradiant phase-transition models, such as the Dicke model [47] and the Tavis-Cummings model [48], can also be mapped onto such a system after a Holstein-Primakoff transformation [49,50]. We concretely compare local and global quantum master equations, with and without rotating-wave approximation for the oscillator-oscillator interaction, to exact results provided by a quantum-Langevin-equation description [42- 45]. We explicitly evaluate the stationary mean occupation number of one oscillator for various nonequilibrium temperature differences. We find that the local master equation generally fails to reproduce the results of the quantum Langevin equation especially for large temperature differences, while the global approach exhibits better agreement. We show that this feature is directly related to the inability of the local description to correctly capture intersystem correlations, which we quantify with the help of the quantum mutual information [9].

II. COUPLED-OSCILLATOR MODEL
We consider a system of two interacting harmonic oscillators with Hamilton operator, where a † j and a j are the usual ladder operators and ω j the respective frequencies. We will examine two different types of intersystem interactions: (i) a position-position interaction, x 1 x 2 , corresponding to κ = λ, and (ii) its rotating-wave version, obtained for κ = 0. Two important points should be stressed: First, the positionposition coupling x 1 x 2 leads to critical behavior above a critical interaction strength [51][52][53], in contrast to the commonly treated Hookian interaction (x 1 − x 2 ) 2 [41][42][43][44][45]. In addition, while the rotation-wave approximation is usually associated with a weak-coupling condition, λ/ω i 1, it has recently been shown that counterrotating terms may be effectively suppressed in modulated systems, even in the ultrastrong regime [54,55]. This opens the possibility to experimentally study critical behavior in strongly interacting rotating-wave models.
The isolated Hamilton operator (1) may be diagonalized exactly for both intersystem interactions, yielding two uncoupled modes with respective energies [29,51,52], These energies display critical behavior at the respective critical couplings λ pp c = √ ω 1 ω 2 /2 and λ rw Above these points, the eigenfrequencies of the Hamilton operator (1) become imaginary or negative. The energy spectrum is thus no longer bounded from below. These critical values are in agreement with those of the Dicke and Tavis-Cummings models [56][57][58]. We next attach each quantum harmonic oscillator to a heat bath with respective temperature T j (Fig. 1). As commonly done, we model these reservoirs by an ensemble of harmonic oscillators [2][3][4][5]. We further assume that the system-bath coupling is weak, so that the rotatingwave approximation is applicable to that coupling (see details in the Supplemental Material [59]). We reemphasize that the interaction between the two harmonic oscillators of the system (1) might be strong.
Most studies of dissipative phase transitions consider complex interacting many-body systems [17][18][19][20][21][22][23][24][25][26][27][28]. The direct comparison between global and local master equation descriptions is thus extremely difficult in these systems. By contrast, the coupled-oscillator model is complicated enough to exhibit steady-state critical behavior and, at the same time, simple enough to allow for (i) a detailed comparison between global and local approaches, and (ii) the evaluation of intersystem correlations.

III. QUANTUM-MASTER-EQUATION DESCRIPTION
In the usual Born-Markov limit, the density operator ρ of the joint quantum system obeys a Lindblad master equation of the form [2-5] (we set = 1 throughout), where the dissipators are given by The coefficients Γ k (A i , A j ), as well as the operators A i , depend on the local or global type of the quantum master equation [29].
In the local approach, each oscillator interacts with its heat bath (labelled by k = 1 or 2) as if it were not coupled to the other oscillator. As a result, the quantum master equation may be derived as usual in the local eigenbasis of one oscillator [2][3][4][5]. The operators A i are here the standard ladder operators, (a 1 , a † 1 , a 2 , a † 2 ), and the dissipators are given by denotes the thermal occupation number [2][3][4][5]. These formulas evidently hold for the two kinds of oscillatoroscillator interaction in Eq. (1).
On the other hand, the global master equation is derived in the global eigenbasis of the combined twooscillator system [29,43]. The diagonalization of the joint Hamilton operator (1) accounts for the indirect subsystem-reservoir and reservoir-reservoir correlations which are generated by their coupling to the system. Such correlations are ignored in the local approach. This is the reason why the local master equation may violate the second law of thermodynamics [29]. The explicit (and lengthy) expressions for the dissipators are summarized for both the position-position and rotating-wave interactions in the Supplemental Material [59]. In this situation, they depend on operators A i that are given by properly rotated ladder operators [59].
In the following, we will solve the four different quantum master equations (global/local forms with/without rotating-wave interaction) by applying a characteristic function method in symplectic space [59] and evaluate the steady-state mean occupation number a † j a j ss = tr(ρ ss a † j a j ), where ρ ss is the stationary density operator.

IV. QUANTUM-LANGEVIN-EQUATION DESCRIPTION
In order to assess its validity, both for equilibrium and nonequilibrium conditions, we shall compare the steadystate properties of the approximate quantum-masterequation treatment to those of the exact quantum-Langevin-equation approach [3]. To this end, we will extend the results obtained for Hookian coupling [42][43][44][45] to the position-position and rotating-wave interactions of Eq. (1). The quantum Langevin equations read [3], where the noisy input operators a j,in , stemming from the interaction with the respective baths, are characterized by the correlation function in Fourier space, The coupled quantum Langevin equations (5) can be solved by matrix inversion in Fourier space [59]. In particular, the mean occupation number is here equal to, (8) is independent of time in the steady-state regime and we will set t = 0 in the following. Steadystate mean occupation numbers may be evaluated exactly (without any approximations) in the quantum-Langevinequation formalism in contrast to the quantum-masterequation approach [42]. Figure 2 presents the steady-state mean occupation number a † 1 a 1 ss of the first oscillator as a function of the reduced interaction strength λ/λ c in the equilibrium (high-temperature) case ∆T = T 2 − T 1 = 0. We observe perfect agreement between the global quantum master equation (blue dots), the quantum Langevin equation (green line) and the equilibrium (Gibbs) state ρ eq = exp(−βH)/Z (yellow line) [59] for all values of λ/λ c , for both the position-position and the rotating-wave (inset) interactions. By contrast, the local quantum master In order to gain deeper insight on the nonequilibrium properties of the different quantum master equations, we next examine the ratio of their steady-state mean occupation numbers and the corresponding quantum-Langevinequation expressions, a † 1 a 1 ss / a † 1 a 1 Langevin , for increasing temperature differences ∆T . In the high-temperature regime (β i ω i 1), Fig. 3a shows that, for the positionposition interaction, the local approach gets worse as the system moves further away from equilibrium and that even the global approach slightly departs from the predictions of the quantum Langevin equation for large ∆T . A similar behavior is seen in Fig. 3b when the two temperatures are low (β i ω i 1). Analogous results are displayed for the rotating-wave interaction in Figs. 4ab: remarkably, the global quantum master equation here always perfectly matches the quantum Langevin equation, for all λ and all ∆T , while the local quantum master equation always fails to describe critical behavior. We also mention that the discrepancy between the various descriptions in general depends on the sign of the nonequilibrium temperature difference ∆T [59].

V. RESULTS
The success/failure of the quantum-master-equation description of dissipative critical phenomena may be understood both physically and mathematically. To first address the physical aspect, we consider the quantum mutual information between the two harmonic oscillators, I(ρ) = S(ρ 1 ) + S(ρ 2 ) − S(ρ), where S(ρ i ) = −tr{ρ i ln ρ i } is the von Neumann entropy and ρ i = tr i ρ are the reduced density operators of the respective harmonic oscillators [9]. The quantum mutual information is a measure of the total (classical and quantum) correlations between two subsystems and has been used broadly to characterize critical transitions [60][61][62][63][64]. Figure 5 shows that the stationary quantum mutual information I(ρ ss ) displays a very similar dependence on the interaction strength λ as the average occupation number a † 1 a 1 ss represented in Fig. 1, both for the positionposition and rotating-wave interoscillator interactions. The shortcomings of the quantum-master-equation approach, especially in its local version, may thus be traced to its inability to correctly capture intersystem correlations close to the critical point. This feature can be confirmed mathematically by looking at the way the respective Lindblad quantum master equations are obtained [29]: the dissipators in the local master equation are indeed derived in the local eigenbasis of each separate harmonic oscillator, while those of the global master equation follow from a diagonalization of the interacting two-oscillator system (the unitary evolution given by the von Neumann term in Eq. (4) describes coupled dynamics in both cases). The global scheme thus better accounts for intersystem correlations than the local one, and should therefore be preferred. Such intersystem correlations are indeed crucial for an accurate description of many-body critical systems, and should not be incorrectly omitted Yet, despite these deficiencies, local quantum master equations have been a tool of choice in numerous studies on dissipative critical behavior [17][18][19][20][21][22][23][24][25][26][27][28].

VI. CONCLUSIONS
We have examined the ability of global and local quantum master equations to accurately describe dissipative critical phenomena using an illustrative system of two interacting, damped harmonic oscillators, with and without rotating-wave interaction. This model provides a transparent, yet generic, example to perform such a study. We have found that while the global master equation reproduces the results of the quantum Langevin equation reasonably well, the local version usually fails to do so, especially in the far-from-equilibrium regime; it generally fails in the case of the rotating-wave interaction. We have related these properties to the inability of the local approach to correctly apprehend oscillator-oscillator correlations that we have quantified with the help of the quantum mutual information. The latter quantity could be easily determined in the present two-oscillator model, in contrast to more complex interacting many-body systems. Our findings show that approximate local quantum master equations in general, and their exact analytical solutions in particular, should be used with caution when studying dissipative critical behavior, and that the more complicated global approach should be favored instead.

ACKNOWLEDGMENTS
We acknowledge financial support from the German Science Foundation (DFG) (Contract No FOR 2724).

APPENDIX A: QUANTUM MASTER EQUATIONS
In the following, we provide details about the dissipators of the four quantum master equations that we consider in our study (global/local forms with/without rotating-wave interaction) as well as their solutions.
The standard dissipators of the local quantum master equation are given below Eq. (4) in the main text. They are derived in the local eigenbasis of each oscillator [2-5] and thus hold for both the position-position and rotating-wave intersystem interactions. By contrast, the global master equations are derived in the global eigenbasis of the combined two-oscillator system obtained by diagonalizing the quadratic Hamilton operator H. The respective dissipators are then computed by expanding the system-bath interaction in this basis. We concretely consider the total system-bath Hamilton operator, with harmonic thermal baths H Bi = j ω ij b † ij b ij and local system-bath couplings H SBi = j κ ij (a i b † ij + h.c.) with coupling constants κ ij (i = 1, 2) [2][3][4][5].
(2) of the main text. The global Lindblad dissipators are then for the quantum oscillator 1 coupled to bath 1 at inverse temperature β 1 with W = S −1 . Here the indexes i, j run over 1-4, corresponding to the elements of (a 1 , a 2 , a † 1 , a † 2 ) and k, l run over all combinations of 1,3. The indexes k, l correspond to the initially chosen local coupling terms in the derivation of the master equation before the diagonalization is applied, which can be ordered either as a i a † i or a † i a i . Thus, there are 32 different terms corresponding to the 16 unique operator orderings A i A j in the dissipators. Expressions for second bath at inverse temperature β 2 are analogous with k, l now combinations of 2,4.
We explicitly solve the linear local and global quantum master equations by computing the first and second moments of ρ in symplectic space [2]. The symmetric characteristic function is defined by is the displacement operator. The (symmetric) moments are then obtained by differentiation [65], where · s is the expectation value of the symmetrized version of the operators a †k i a l j . The evolution of the characteristic function is derived from the master equation together with the identities, with α i = x i + ip i and d/dα i = (d/dx i − id/dp i )/2 using the Gaussian ansatz χ(x 1 , p 1 , x 2 , p 2 ) = exp(i P ȳ − P Tσ P /2) with P = (x 1 , p 1 , x 2 , p 2 ) T and ȳ = (ȳ 1 ,z 1 ,ȳ 2 ,z 2 ) T . Since the Hamiltonian is purely of quadratic order, the steady state values for the first moments always vanishȳ i = 0 =z i and the system is completely described by the second moments.

APPENDIX B: QUANTUM MUTUAL INFORMATION
The quantum mutual information for a Gaussian system can be calculated from the covariance matrix as p 1 , x 2 , p 2 ). In this form, the submatrices of interest are σ = ((α, γ), (γ T , β)).
The Gibbs state expectation values may in addition be evaluated by the diagonalization of the Hamiltonian given above. In general, any quadratic expectation value in the a i algebraic space may be calculated via A i A j = k S ik S j C k C with A i = (a 1 , a 2 , a † 1 , a † 2 ) and C i = (c 1 , c 2 , c † 1 , c † 2 ). The C i C j are then given by the uncoupled oscillators with eigenfrequencies, Eq. (3) of the main text, and corresponding temperature T .

APPENDIX D: NONEQUILIBRIUM STEADY STATES
The deviation of the local (and, to a lesser extent, global) quantum master equations from the quantum Langevin equation does not only depend on the magnitude of the temperature difference ∆T but also on its sign (Figs. 6ab). We first note that the local master equation leads to larger (smaller) mean occupation numbers for weak (strong) coupling, both for the positionposition and the rotating-wave interactions. This increase of the mean occupation number at small coupling is caused by the frequency difference between the oscillators (ω 1 > ω 2 ), which leads to a relatively larger occupation number in the second oscillator (modulated by the temperature difference), whereas its decrease is induced by strong-coupling effects. In addition, the global master equation completely matches the Langevin equation, for all ∆T for the rotating-wave interaction, while this is not the case for the position-position interaction: the mean occupation number is larger (smaller) than that the quantum Langevin for ∆T < 0 (∆T > 0).