Research ArticlePHYSICS

# General description for nonequilibrium steady states in periodically driven dissipative quantum systems

See allHide authors and affiliations

Vol. 6, no. 27, eabb4019

## Abstract

Laser technology has developed and accelerated photo-induced nonequilibrium physics, from both the scientific and engineering viewpoints. Floquet engineering, i.e., controlling material properties and functionalities by time-periodic drives, is at the forefront of quantum physics of light-matter interaction. However, it is limited to ideal dissipationless systems. Extending Floquet engineering to various materials requires understanding of the quantum states emerging in a balance of the periodic drive and energy dissipation. Here, we derive a general description for nonequilibrium steady states (NESSs) in periodically driven dissipative systems by focusing on systems under high-frequency drive and time-independent Lindblad-type dissipation. Our formula correctly describes the time average, fluctuation, and symmetry properties of the NESS, and can be computed efficiently in numerical calculations. This approach will play fundamental roles in Floquet engineering in a broad class of dissipative quantum systems from atoms and molecules to mesoscopic systems, and condensed matter.

## INTRODUCTION

State-of-the-art laser technology has opened new research fields in physics: the Floquet science and engineering (13). The main focus of these fields is the nonequilibrium states driven periodically by external fields, e.g., intense laser fields. Physical properties of the nonequilibrium states are mainly understood by the so-called effective Hamiltonian, which reflects the periodic driving, according to the Floquet theorem (4) and the ensuing theoretical developments (58). Conversely, designing a suitable driving protocol, one can engineer the effective Hamiltonian, which enables us to have desirable properties and functionalities of physical systems. Various exotic states and useful manipulation of matter have been theoretically proposed, and some of them have been experimentally realized: Floquet topological states (9) in solids (10), ultracold atomic gases (11), and in photonic wave guides (12), Floquet time crystals (13) in nitrogen-vacancy (NV) centers (14) and trapped ions (15), and control of quantum magnets (16) and their interactions (17).

However, these Floquet-theoretical predictions based on the effective Hamiltonian are quantitatively reliable only in ultraclean materials or well-designed artificial systems, where dissipation is negligible. For the Floquet science and engineering in a variety of materials, it is indispensable to understand the nonequilibrium steady state (NESS), which emerges in a balance of the energy injection by the periodic driving and the energy dissipation (1821). For individual systems, by considering specific sources of dissipation, i.e., system-bath couplings, one can calculate physical quantities in the NESS and predict interesting phenomena such as Floquet topological insulators (22, 23), periodic thermodynamics (24), dynamical localization (25), and generalized Bose-Einstein condensation (26). In this research direction, the Floquet-Green function approach has developed and enables us to calculate various physical effects dependent on the type of the system-bath coupling (2729). Another research direction, which we address here, is to seek for a universal characterization for the NESS. We could imagine that there exists a simple and general expression for the NESS when the dissipation is weak and featureless. An attempt is to conjecture that the NESS is generally described by the Floquet-Gibbs state (FGS) (21, 30, 31), i.e., the Gibbs state with the effective Hamiltonian, but the conditions for the FGS being realized have shown quite restrictive (30). Hence, despite its importance, the general formula for the NESS has been still an elusive problem.

Here, in exchange for restricting ourselves to the high-frequency drivings, we deal with generic systems and driving protocols, obtaining simple and general formulas for the NESS (Eqs. 9 to 11 below). We derive these formulas by applying the high-frequency expansion technique, which has been recently developed (5, 6, 3234), to the Lindblad equation with periodic Hamiltonians. As exemplified in an effective model for the NV center in diamonds (35), our formulas correctly describe both the time average and fluctuation of the NESS at the leading order of ω−1 (ω denotes the driving frequency). These formulas also capture nontrivial behaviors of physical quantities due to the dynamical-symmetry breaking that cannot be described by the effective Hamiltonian or the FGS and will thereby play critical roles in the Floquet science and engineering in dissipative quantum systems.

## RESULTS

### Floquet-Lindblad equation with detailed balance condition

We begin by considering a quantum system defined on an N-dimensional Hilbert space. This system can be single body or many body as long as it satisfies the requirements that will be described below. We let H0 denote the time-independent Hamiltonian, which describes our system in the absence of driving. The eigenenergies and eigenstates of H0 are denoted by {Ei}i=1N and {Ei}i=1N, respectively. For simplicity, we assume that the eigenenergies are not degenerate and E1 < E2 < ⋯ < EN (the generalization to degenerate H0 is formulated in the Supplementary Materials). The effect of the driving is represented by a time-dependent part Hext(t) of the total HamiltonianH(t)=H0+Hext(t)(1)We assume that the driving term is periodic with period T: Hext(t + T) = Hext(t) and hence H(t + T) = H(t). Without loss of generality, the decomposition (Eq. 1) is defined so that the time average of Hext(t) vanishes, 0Tdt Hext(t)=0. Thus, the Fourier series of Hext(t) can be written asHext(t)=m0Hmeimωt(2)

To study driven dissipative systems, we consider the density operator ρ(t), whose dynamics is described by the Floquet-Lindblad equation (3639) (we set ℏ = 1 throughout this paper)dρ(t)dt=Ltρ(t)=i[H(t),ρ(t)]+D[ρ(t)](3)D[ρ(t)]i,jΓij[Lijρ(t)Lij12{LijLij,ρ(t)}](4)Here, Lij ≡ ∣Ei〉〈Ej∣ is the time-independent Lindblad operator describing the transition from the j-th to the i-th eigenstates of the undriven Hamiltonian H0. Each Lindblad operator Lij represents a decay (excitation) process for i < j (i > j). The real number Γij ( ≥ 0) denotes the rate for the corresponding process, and we set Γii = 0 for each i. The transition rates Γij must be small enough for the Floquet-Lindblad equation to be valid (see Discussion). Note that Eq. 3 is trace-preserving d tr[ρ(t)]/dt = 0, and thus, we use the normalization tr[ρ(t)] = 1. We assume that the transition rates Γij satisfy the detailed balance conditionΓijeβEj=ΓjieβEi (for ij)(5)where β is the inverse temperature of the bath coupled to the system. We also assume that the matrix Γij is a nonnegative irreducible matrix (40). These assumptions ensure that, without driving, the system goes, irrespective of the initial state, to the thermal equilibrium state, or the canonical ensemble ρcan = e−βH0/Z of H0 with Z = tr(e−βH0). We note that the Lindblad operators Lij may depend on the driving in general if we consider more microscopic theories of dissipation (21). However, we neglect this dependence in this work for simplicity.

### General formulas for NESS

The key idea to obtain the NESS is the high-frequency expansion for the Floquet-Lindblad equation (34). Among several formulations, we adopt the van Vleck perturbation theory (5, 6), which leads to the following propagation for ρ(t) (see Materials and Methods): ρ(t)=eG(t)e(tt)LeffeG(t)ρ(t). The time-independent part ℒeff is represented by the effective HamiltonianLeff(ρ)=i [Heff,ρ]+D(ρ)+O(ω2)(6)with Heff = H0 + ω−1n > 0[Hn, Hn]/n + O−2). The time-dependent part eG(t) is the so-called micromotion operator periodic in time G(t+T)=G(t) and is given by G(t)(ρ)=ω1m0[Hm,ρ]eimωt+O(ω2). Without loss of generality, we suppose the initial time to be t′ = 0, havingρ(t)=eG(t)etLeffeG(0)ρ(0)(7)with ρ(0) being our initial state.

To obtain the asymptotic behavior of ρ(t), we focus on the first part ρ(t)=etLeffeG(0)ρ(0). Note that this is the solution of the time-independent Lindblad equation dρ′(t)/dt = ℒeff ρ′(t) from the initial state eG(0)ρ(0). Under our assumptions on D, ρ′(t) approaches, irrespective of the initial state, the unique state ρ characterized by Leffρ=0. Thus, we come to the first main result, obtaining the asymptotic behaviorρ(t)ρness(t)=eG(t)ρast(8)Because G(t)=G(t+T), this NESS is also periodic in time. Focusing on the leading-order contribution, we have a simple explicit formula for ρness(t)ρness(t)=ρcan+σMM(t)+σFE+O(ω2)(9)in which both σMM(t) and σFE are O−1), and we call σMM(t) and σFE the micromotion and Floquet-engineering parts, respectively. Equation 9 is our second main result, which we prove in Materials and Methods. Its generalization in the absence of the detailed balance condition is outlined in Discussion. Note that tr[ρness(t)] = 1 is satisfied, at least, up to this order because both σMM(t) and σFE are traceless, as will be evident below.

The micromotion part σMM(t) is defined byσMM(t)=1ωm0eimωtm[Hm,ρcan](10)We have named it after the following two properties of σMM(t). First, this part is periodic in time σMM(t + T) = σMM(t) and contributes to oscillations of physical observables. Second, σMM(t) does not contribute to the time averages of physical observables for one period of oscillation. For an observable A, we have t1t1+Tdt tr[σMM(t)A]/T=0.

The Floquet-engineering part σFE in Eq. 9 is independent of time and given byEkσFEEl=EkΔHeffEl(EkEl)iγkl(pcan(k)pcan(l)) (for kl)(11)and 〈Ek∣σFEEk〉 = 0 for all k, where ΔHeffHeffH0 = O−1), pcan(k)=eβEk/Z is the Boltzmann weight, and γkli(Γik+Γil)/2 represents the symmetric transition-rate matrix (see the Supplementary Materials for the generalization to degenerate H0). The reason why we call σFE the Floquet-engineering part is that it describes how the effective Hamiltonian changes physical observables from their values in thermal equilibrium. In contrast to the micromotion part, the Floquet-engineering part contributes to the time-averaged quantities. As we will show below, Eq. 11 is regular in the weak dissipation limit γij → 0, where ρness(t) becomes independent of γij and coincides with the canonical Floquet steady state (CFSS) that we define in Eq. 14 below.

Equation 11 serves as the foundation for the Floquet engineering in dissipative quantum systems. Let us imagine, for example, that an observable A has zero expectation value at thermal equilibrium, tr(ρcanA) = 0, but nonzero value for the NESS, tr(σFEA) ≠ 0. This situation means that one can implement an appropriate periodic driving Hext(t) and hence ΔHeff, thereby activating the observable A. Upon this engineering, Eq. 11 tells us how much activation is possible for observables of interest. We will see some examples below.

In addition to their generality, our formulas (Eqs. 9 to 11) are extremely efficient in practical calculations of the NESS. In the straightforward numerical calculation, one integrates the Floquet-Lindblad equation (Eq. 3) with a sufficiently short time step until the system reaches the NESS. In contrast, our formulas enable us to evaluate the NESS at an arbitrary time t without numerical integration once we have the energy eigenstates {∣Ek〉} of the time-independent Hamiltonian H0. This difference of efficiency becomes more evident when the Hilbert-space dimension N is large. In the straightforward calculation, the density matrix ρ(t) is commonly treated as an N2-dimensional vector and the superoperator ℒt as an N2 × N2 matrix. Thus, the computational complexity for each time step is O(N4) in general. On the other hand, the complexity for our formula is one order smaller and given by O(N3), which derives from the exact diagonalization of H0. Thus, our formulas enable us to evaluate the NESS for larger Hilbert-space dimensions that occur, for example, in quantum many-body systems. For special cases where ℒt is a sparse matrix and has only O(N0) nonzero elements, the computational complexity for one time step is O(N2). Nevertheless, even for these cases, we need many time steps typically larger than N in obtaining accurate results, and hence, our formulas require less computational complexity.

### Numerical verification with the NV center in diamonds

By taking a single spin with S = 1, we demonstrate how our formula (Eq. 9) works in the quantum dynamics described by the Floquet-Lindblad equation (Eq. 3). We consider an effective Hamiltonian for the NV center in diamonds (35)HNV(t)=BsSz+NzSz2+Nxy(Sx2Sy2)+Hextcirc(t)(12)where Bs is the static Zeeman field, Nz and Nxy are the coupling constants of magnetic anisotropic terms, and Hextcirc(t)Bd(Sxcosωt+Sysinωt) represents the coupling to the circularly polarized ac magnetic field. We note that the energy eigenstates of the time-independent part of Eq. 12 are analytically obtained, and our formula can be computed almost analytically in this model. Because NzNxy in the NV centers (35), we set Nz=1 and Nxy=0.05 in our analysis. The ac-field irradiation leads to spin dynamics, as schematically illustrated in Fig. 1A.

Typical time evolutions A(t) ≡ tr[ρ(t)A] are shown for two representative observables A = Sz and {Sx, Sy} ( = SxSy + SySx) in Fig. 1 (B and C). In this calculation, we take the thermal state for the time-independent part of HNV at t = 0 and let the state evolve according to the Floquet-Lindblad equation (Eq. 3) with H(t) being HNV(t). The static Zeeman field is Bs = 0.3, and the driving parameters are Bd = 0.1 and ω = 1. As for the Lindblad operators, we take Γij according to the heatbath method as Γij = γeβEi/(eβEi + eβEj) for ij, with rate constant γ = 0.2 and inverse temperature β = 3, and Γii = 0 for all i’s. Figure 1 (B and C) shows that after a sufficiently long time t ≫ γ−1, the system reaches the NESS, in which the observables oscillate with period T = 2π/ω. In particular, the observable A = {Sx, Sy} is initially zero for a symmetry reason (e.g., Sy → − Sy) but becomes nonzero on time average. Namely, this observable is engineered by the periodic drive Hextcirc(t).

To test our formula (Eq. 9) quantitatively, we first focus on the one-cycle average A¯(ω)=T1tt+Tds tr[ρ(s)A] for t ≫ γ−1. For this test, we also define A¯HF(ω)=T1tt+Tds tr[ρHF(s)A], where ρHF(s) is the density matrix calculated from Eqs. 9 to 11 based on the high-frequency expansion. In Fig. 1 (D and E), we compare the one-cycle averages A¯(ω) and A¯HF(ω) [recall that the micromotion part σMM(t) does not contribute to the one-cycle averages]. At high-frequency ω ≳ 10, the difference ΔAHF(ω)=A¯HF(ω)A¯(ω) decreases quite well. In Fig. 2 (A and B), we plot ΔAHF(ω) against ω for A = Sz and {Sx, Sy}, respectively. We stress that the difference ΔAHF(ω) decreases more rapidly than O−1) [O−2) for Sz and O−3) for {Sx, Sy}]. This means that Eqs. 9 and 11 perfectly describe the actual NESS at the level of O−1). As shown in the Supplementary Materials, ΔAHF(ω) = O−2) holds true not only for these two observables but also for all the other observables. Therefore, we have verified Eq. 9 apart from the micromotion part.

For the complementary test of our formula (Eq. 9), we consider the one-cycle standard deviation (SD) ΣA(ω)=[T1tt+Tds{tr[ρ(s)A]A¯(ω)}2]1/2, which quantifies the micromotion amplitude. Similarly to A¯HF(ω), we also introduce the one-cycle SD based on our high-frequency formulas, ΣAHF(ω)=[T1tt+Tds{tr[ρHF(s)A]A¯HF(ω)}2]1/2. These quantities are suitable for testing Eq. 10 because they include the contribution only from the micromotion part σMM(t). Because ΣA(ω) is an O−1) quantity in general, the accuracy of our formula is verified if the difference ΔΣAHF(ω)=ΣAHF(ω)ΣA(ω) is O−2). This criterion is satisfied, as shown in Fig. 2 (C and D) for A = Sz and {Sx, Sy}, respectively. We remark that our formula leads to ΣA(ω) = 0 at O−1) for these observables, which can be analytically shown by noticing H±1S± = Sx ± iSy. Thus, the plotted data correspond to ΣA(ω) itself for the actual dynamics, and ΔΣA(ω) could be reduced by dealing with the higher-order terms in Eq. 9. In any case, the fact that ΔΣA(ω) is O−2) justifies Eqs. 9 and 10.

### Comparison with the FGS

Let us make comparisons with the FGS, which has been a candidate for the NESS description in periodically driven dissipative quantum systems (21, 30, 31). To define the FGS, we introduce the Floquet state ∣ui(t)〉 and its quasienergy ϵi. According to the Floquet theorem, the time-dependent Schrödinger equation iddtψ(t)=H(t)ψ(t) has the independent solutions ∣ψi(t)〉 = eiϵitui(t)〉 (i = 1,2, …, N), with periodicity ∣ui(t + T)〉 = ∣ui(t)〉 and 〈ui(t)∣uj(t)〉 = δij. In terms of the Floquet states, the FGS is defined byρFG(t)=1ZFGieβϵiui(t)ui(t)(13)where ZFG=ieβϵi. To obtain the Floquet states and quasienergies in practice, the common method, which we use here, is to calculate the one-cycle unitary evolution U(T)=Texp[i0Tds H(s)], where Texp denotes the time-ordered exponential, by numerical integrations of the time-dependent Schrödinger equation. The eigenvectors and eigenvalues of U(T) correspond to ∣ui(0)〉 and eiϵiT, which give us ϵi and ∣ui(t)〉. Note that the Floquet states and quasienergies thus obtained are exact and involve all-order contributions in 1/ω.

Quantitative comparisons between the actual dynamics and the FGS are shown in Fig. 2. For the FGS, we define the one-cycle average A¯FG(ω)=T1tt+Tds tr[ρFG(s)A] and the SD ΣAFG(ω)=[T1tt+Tds{tr[ρFG(s)A]A¯FG(ω)}2]1/2. In Fig. 2 (A and D), we plot the difference ΔAFG=A¯FG(ω)A¯(ω) for the two representative observables A = Sz and {Sx, Sy}. Remarkably, the difference is O−1), meaning that the FGS cannot reproduce the leading-order contribution of the one-cycle average in contrast to our Eq. 9. As for the micromotion amplitude, or the one-cycle SD ΣA(ω), Fig. 2 (C and D) shows that the FGS reproduces the actual values better than our Eqs. 9 and 10. This is, in part, because the FGS involve all-order contributions in 1/ω, while our formulas are the leading-order approximation. As stated above, our formula could be improved order by order when we start over from our first main result (Eq. 8).

A weak point of the FGS is highlighted in Fig. 1E, in which the FGS gives {Sx,Sy}¯(ω)=0 for any ω, while it is not true in the actual dynamics. This is due to an antiunitary dynamical symmetry constraining the Floquet states and, hence, the FGS. We take an antiunitary operator V: VSyV = − Sy and VSαV = Sα (α = x and z). Then, we notice the dynamical symmetry VHNV(Tt)V = HNV(t), which implies that u˜i(t)Vui(Tt) is also a Floquet state with quasienergy ϵi. Assuming that quasienergies are not degenerate as in our examples, we have that u˜i(t) and ∣ui(t)〉 are equivalent up to an overall phase shift. Because of VAV = − A with A = {Sx, Sy}, the one-cycle averages of A calculated for u˜i(t) and ∣ui(t)〉 differ by their signs, meaning that the one-cycle average vanishes. Note that similar arguments apply to other observables, satisfying VAV = − A.

We remark that the dissipation can break such an antiunitary dynamical symmetry, and this is the origin of the nonzero one-cycle average of {Sx,Sy}¯(ω). We can show that this average vanishes by taking the limit γij → 0 in Eq. 11. In other words, the NESS in dissipative systems shows richer properties inferred only from the effective Hamiltonian itself. Our Eq. 9 well describes these properties unlike the FGS (Eq. 13), which incorporates no information about the dissipation, or the Lindblad operators.

One might be interested in an approximate description of the NESS independent of the details of γij for weak dissipation and expect that the FGS serves such a description. This is not true, at least within our formulation of periodically driven dissipative systems described by Eqs. 3 and 4. Instead, the actual NESS coincides with yet another state, which we name the canonical Floquet steady state (CFSS), defined byρCFSS(t)=1ZieβEiui(t)ui(t)(14)Note that the CFSS is similar to the FGS (Eq. 13) and obtained by replacing the quasienergy ϵi with the real energy Ei. One can analytically show that the CFSS at high frequency coincides with our formulas (Eqs. 9 to 11) in the limit of γij → 0 (see the Supplementary Materials for details).

## DISCUSSION

We have derived and verified the simple and general formulas (Eqs. 9 to 11) describing the NESS in dissipative quantum systems under high-frequency periodic drive. We have also exemplified the breaking of an antiunitary dynamical symmetry and the possibility of the Floquet engineering in NV centers in diamonds. Being quite general, our formulas should apply to various quantum systems such as atoms and molecules, trapped ions, condensed matters, and so on, and play the fundamental role in understanding and engineering unusual nonequilibrium states in these systems.

The parameter region in which our formulas are valid is depicted in Fig. 3. On the basis of the high-frequency expansion, our formulas are valid when the driving frequency ω (more precisely, the photon energy ℏω) is greater than the energy scales of the system and the system-drive coupling. We note, however, that our formulas hold true for any strength of dissipation Γij (or γij) within the Floquet-Lindblad equation (Eq. 3). As we have shown, our formulas reduce to the CFSS rather than the FGS when the dissipation strength is smaller than the energy gap, i.e., the nonzero minimum difference between eigenenergies (our formulas are generalized for the degenerate Hamiltonian in the Supplementary Materials).

One should note that the Floquet-Lindblad equation (Eq. 3) becomes invalid when Γij is too large. The Lindblad-type dissipation is derived from several approximations such as the Born-Markov approximation (39). These approximations require the condition that the relaxation time ~1/Γij is longer than the time scale of the system’s dynamics and the bath correlation time. Namely, the dissipation rate Γij should be smaller than the other relevant energy scales. Note that the high-frequency driving does not break this condition, while lower frequencies may be problematic.

We remark a further generalization of our results Eqs. 9 to 11. Although we have assumed the detailed balance condition (Eq. 5), this condition can be removed as long as the transition-rate matrix Γij is irreducible. In this case, the solution of the Lindblad equation without driving is not the canonical ensemble ρcan but another state ρ˜ characterized by i[H0,ρ˜]+D(ρ˜)=0. Correspondingly, our results for the NESS (Eqs. 9 to 11) hold true with the following replacements: ρcanρ˜ and pcan(k)p˜(k)=Ekρ˜Ek. Thus, our formulas apply to any periodically driven dissipative systems as long as the dissipation is of Lindblad type and irreducible. Therefore, our formulas are useful for a broad class of systems in exploring generic features of the NESS and in estimating Floquet-engineered physical quantities.

It remains an open question to find a simple and general formula for the NESS at lower frequency. The applicability of the CFSS in many-body systems is also a nontrivial issue because the energy gap can be very small in those systems. Addressing these issues will lead us to the complete understanding of the NESS in dissipative Floquet systems.

## MATERIALS AND METHODS

### High-frequency expansion for the Floquet-Lindblad equation

The Floquet-Lindblad equation that we discuss in this work is symbolically represented as ∂tρ(t) = ℒ(t)ρ(t), where the time-dependent Liouvillian ℒ(t) is defined byL(t)ρ=i[H(t),ρ]+D(ρ)(15)We introduce the Fourier series for the Liouvillian as ℒ(t) = Σmmeimtω. Because the Lindblad operators Lij are time independent in this work, each Fourier component is given as followsL0ρ=i[H0,ρ]+D(ρ) (16)Lmρ=i[Hm,ρ] (form0)(17)The formal solution of the Floquet-Lindblad equation is obtained as ρ(t)=V(t,t)ρ(t), with the propagator V(t,t)=Texp[ttL(s)ds], where Texp denotes the time-ordered exponential. The determining equations for V are tV(t,t)=L(t)V(t,t) and V(t,t)=1.

The high-frequency expansion in terms of the van Vleck approach makes the following ansatzV(t,t)=eG(t)e(tt)LeffeG(t)(18)where G(t) is periodic in time and ℒeff is time independent. This ansatz satisfies one determining equation V(t,t)=1 for any choices of G(t) and ℒeff, and what determines these two is tV(t,t)=L(t)V(t,t). As we will see below, this equation only determines the derivative of G(t), and thus, we further impose 0TG(t)dt=0 to fix the constant of integration.

To obtain the determining equations for G(t) and ℒeff, we substitute Eq. 18 into tV(t,t)=L(t)V(t,t), having t(eG(t))+eG(t)Leff=L(t)eG(t). With some algebra (see the Supplementary Materials), we rewrite this equation, obtainingtG(t)=k=0Bkk!(adG)k[L(t)+(1)k+1Leff](19)where adG is defined by adGρ=[G(t),ρ], and Bk denotes the Bernoulli number (B0 = 1, B1 = − 1/2, B2 = 1/6, …).

We determine G(t) and ℒeff from Eq. 19 by the series expansionsG(t)=k=1G(k)(t);Leff=k=1Leff(k)(20)We substitute these expansions into Eq. 19 and find the order-by-order solutions, where we assign an order 1 for ℒ(t) and k for G(k)(t) and Leff(k) [see (6) for the case of unitary dynamics]. With some straightforward calculations (see the Supplementary Materials), we obtain the following low-order solutionsLeff(1)=L0(21)G(1)(t)=iωm0eimωtmLm(22)Leff(2)=iωm>0[Lm,Lm]m(23)Note that higher-order components are O−2), and G(1)(t) is equivalent to that used in the main text. We also note that Leff=Leff(1)+Leff(2) up to this order can be represented by the effective HamiltonianHeff=H0+1ωm>0[Hm,Hm]m+O(ω2)(24)By invoking the Jacobi identity [A, [B, C]] + [B, [C, A]] + [C, [A, B]] = 0, we obtain Eq. 6 in the main text.

### Solution for the NESS

Here, we derive Eqs. 9 to 11 from Eq. 8 in the main text. For this purpose, we solve Leffρ=0 for ρ at the leading order of ω−1. It is convenient to work in the energy eigenbasis and separate the diagonal and off-diagonal partsρ=ρ(d)+ρ(od)(25)ρ(d)=kρkkEkEk(26)ρ(od)=k,l(kl)ρklEkEl(27)

First, we consider the off-diagonal elements of both sides of Leffρ=0, having, for kl[i(EkEl)γkl]ρkli(ρllρkk)EkΔHeffElEk[ΔHeff,ρ(od)]El=0(28)where γkli(Γik+Γil)/2 and ΔHeffHeffH0 = O−1). Equation 28 is transformed asρkl=EkΔHeffEl(EkEl)iγkl(ρkkρll)Ek[ΔHeff,ρ(od)]El(EkEl)iγkl(29)Note that the denominators (EkEl) − iγkl do not vanish because γkl > 0 is ensured by the nonnegativity and irreducibility of Γij. Now, as a working hypothesis, we suppose that the diagonal elements ρkk are O0) as verified later. Then, the first term on the right-hand side of Eq. 29 is O−1) because ΔHeff = O−1). Notice that the second term depends only on the off-diagonal elements ρkl, and Eq. 29 can be solved recursively. This yields the ω−1 expansion for the off-diagonal elements ρkl, whose leading order contribution is given byρkl=EkΔHeffEl(EkEl)iγkl(ρkkρll)+O(ω2)(30)

Next, we consider the diagonal elements of both sides of Leffρ=0, havingEkLeffρEk=iEk[ΔHeff,ρ(od)]El+Σl(ΓklρllΓlkρkk)=0(31)We note Ek[ΔHeff,ρ(od)]El=O(ω2) because both ΔHeff and ρ(od) are O−1). Thus, the diagonal elements ρkk are determined up to O−1) by the equationΣl(ΓklρllΓlkρkk)=0(32)According to the irreducibility and the detailed balance condition of Γkl, we have the unique solution asρkk=pcan(k)=eβEkZ(33)where the error is O−2). This result meansρ(d)=ρcan+O(ω2)(34)Because these diagonal elements ρkk are O0), the working hypothesis introduced above has been verified. By substituting Eq. 33 into Eq. 30, we have the leading-order expression for the off-diagonal elementsρkl=EkΔHeffEl(EkEl)iγkl(pcan(k)pcan(l))+O(ω2)=EkσFEEl+O(ω2)(35)which impliesρ(od)=σFE+O(ω2)(36)We remark that tr(σFE) = 0 because each of the diagonal elements of σFE vanishes.

Now that we have obtained the leading-order expression for ρ=ρ(d)+ρ(od), let us calculate the time-dependent density matrix by ρ(t)=eG(t)ρ. By noticing that ρ(d) is O0) and ρ(od) is O−1) and using the Taylor expansion eG(t)=1+G(t)+O(ω2), we obtainρ(t)=ρcan+σMM(t)+σFE+O(ω2)(37)where we have definedσMM(t)=G(t)[ρ(d)]=1ωm0eimωtm[Hm,ρcan](38)Thus, we have derived the Eqs. 9 to 11 here. We remark tr[σMM(t)] = 0, which follows from the cyclic property of trace, and hence tr[ρ(t)] = 1, at least up to this order.