Response theory for nonequilibrium steady-states of open quantum systems

We introduce a response theory for open quantum systems within nonequilibrium steady-states subject to a Hamiltonian perturbation. Working in the weak system-bath coupling regime, our results are derived within the Lindblad-Gorini-Kossakowski-Sudarshan formalism. We find that the response of the system to a small perturbation is not simply related to a correlation function within the system, unlike traditional linear response theory in closed systems or expectations from the fluctuation-dissipation theorem. In limiting cases, when the perturbation is small relative to the coupling to the surroundings or when it does not lead to a change of the eigenstructure of the system, a perturbative expansion exists where the response function is related to a sum of a system correlation functions and additional forces induced by the surroundings. Away from these limiting regimes however, the secular approximation results in a singular response that cannot be captured within the traditional approach, but can be described by reverting to a microscopic Hamiltonian description. These findings are illustrated by explicit calculations in coupled qubits and anharmonic oscillators in contact with bosonic baths at different temperatures.


I. INTRODUCTION
Response theory provides a means of characterizing the dynamical behavior of systems subject to small perturbations. For physical systems obeying detailed balance, the fluctuation-dissipation theorem relates a system's response to the underlying microscopic correlations [1]. Such fluctuation-response relations form the basis of molecular spectroscopy [2][3][4] and are encoded in properties like the conductivities of transport devices and sensitivities of sensors [5][6][7]. For open quantum systems, connections between response and molecular degrees of freedom are generally unknown in instances away from thermal equilibrium. Given the ubiquitous presence of open quantum systems, it is desirable to uncover such relations. Here we work within the Lindblad-Gorini-Kossakowski-Sudarshan (LGKS) description of an open quantum system and develop a response theory for dynamic perturbations [8][9][10]. We find that it is not typically possible to write the response function as a simple correlation function within the system, due to the weak coupling approximation and the reduced description of the LGKS formalism. Nevertheless, we are able to find practical general routes to response functions and identify special cases where Dyson-like expansions of the propagator still exist, resulting in generalized fluctuationdissipation relations.
In the linear response regime for systems obeying detailed balance, response theory has revealed the * amikamlevy@gmail.com † eran.rabani@berkeley.edu ‡ dlimmer@berkeley.edu fluctuation-dissipation theorems that establish a relation between a system's response to a small perturbation and a time correlation function of certain observables evaluated at equilibrium. For a closed quantum system that follows a unitary dynamics, the linear response function of an observable A to a weak perturbationδ(t)V (t) is given by the well known Kubo formula [11] where V is the perturbation operator added to the original Hamiltonian andδ(t) is a small time dependent parameter, and is Planck's constant. The brackets, . . . , denote average over the stationary distribution of the original Hamiltonian. Kubo's result has a classical analogue uncovered originally by Callen and Welton, [12] and has been extended to stochastic nonequilibrium steady states (NESS's) [13][14][15][16] and nonlinear response [17][18][19][20]. Recent work has also established relationships between these results and the entropy production within stochastic thermodynamics [21]. To extend these results to quantum systems that continuously interact with their environment, either the environment needs to be consider explicitly or a suitable reduced description found. While the former are well suited to accounting for instances of strong system-bath coupling [22][23][24][25][26], describing nonequilibrium steady-states in such formalisms are cumbersome, as the composite system plus environment represents a closed system incapable of dissipating energy. If the couplings between the system and environment can be assumed to be weak, a reduced description using the LGKS master equation is appropriate, and a linear response theory can be formulated in the reduced Liouville space [27][28][29][30][31][32][33] studies along these lines have thus far been largely phenomenological.
In this work, we introduce a response theory for open quantum systems within the LGKS formalism that is consistent with an underlying Hamiltonian description of the full system, sub-system plus environment. In particular, we study the response of an open quantum system at equilibrium or NESS to a small Hamiltonian perturbation. The initial NESS of the system corresponds to a steady-state of a system interacting with multiple baths. The system, which remains in contact with the baths at all times, is then subject to a stationary Hermitian perturbation that drives it away from its initial steady-state, as illustrated in Fig. 1.
We find that even a small Hamiltonian perturbation with respect to the system Hamiltonian may result in nonlinear response, and that unlike the case for closed quantum systems or classical systems, the corresponding response function is not simply related to a correlation function of the system. This occurs when the perturbation is considered small with respect to the system Hamiltonian, but not necessarily small with respect to the system-bath Hamiltonian. In such cases, additional timescales emerge within the system whose implications to the Markovian and secular approximations in the LGKS formalism must be considered.
In cases where the perturbation is small with respect to the system-bath Hamiltonian, or when the perturbation influences only the eigenvectors of the system Hamiltonian, a linear response relation can be established. However, we find that the linear response function is subject to corrections over the closed system result that can be traced to noncommutative effects between the perturbation and the system-bath coupling Hamiltonians. Only with the inclusion of these additional terms is the response thermodynamically consistent.
The observation that the perturbation affects the system-bath couplings and leads to corrections in the response function is closely linked to the local and global approaches to deriving the LGKS master equation [34]. In the local approach, the dynamics of a quantum system composed of several subsystems, in which each is coupled to its own surroundings, is described by a local dissipator that ignores the couplings between the subsystems. The global approach takes the inter-couplings between the subsystems into consideration when deriving the dissipative part. The validity of the different approaches has been studied extensively in recent years and is an ongoing inquiry [34][35][36][37][38][39][40][41][42][43][44][45][46]. This study offers a systematic method to calculate the response functions of open quantum systems, and, in doing so, bridges the two approaches. In particular, using a perturbative treatment, we extend the validity of the local approach without the need of a full global treatment, which may become complicated and, in most cases, not analytically feasible.
The manuscript is laid out in the following way. First in Sec. II, we briefly review the assumptions invoked in the LGKS framework, pay particular attention to the coarse-graining required to satisfy the Markovian approximation. Then in Sec. III, we develop a response theory for impulsive system perturbations. This is specialized to instances where only the eigenvectors change in Sec. III A, with an explicit example of such a theory for two anharmonic oscillators are coupled to two heat baths considered in Sec. III B. A response theory specialized to instances where only the eigenvalues change is shown in Sec. III C, with an explicit example of such a theory for two qubits coupled to two heat baths considered in Sec. III D. The response function of thermodynamic functions that are state dependent is discussed in Sec. IV. Finally, in Sec. V we summarize and discuss the main results.

II. FRAMEWORK FOR THE LGKS MASTER EQUATION
We consider a system-bath model with total Hamiltonian H of the form where the quantum system of interest is described by the unperturbed Hamiltonian, in its eigenbasis. While the baths may exist in the thermodynamic limit, we consider explicitly system Hilbert spaces that are discrete. We focus our analysis on step perturbationsδ(t) = δθ(t), and δ is a unitless parameter that sets the scale, θ(t) is the Heaviside function turned on at t = 0, and V an arbitrary operator in the system Hilbert space. The system is weakly coupled to two or more baths held in distinct thermal states that can generate flows of mass and energy through the system described by where r † α,j r α,j is the number operator for the jth mode of the αth bath, which is assumed have a continuous spectrum and obey bosonic statistics. We will assume that the interaction with a bath is bilinear such that the interaction Hamiltonian is where for the αth bath an operator in the system, S α , is coupled to a bath operator, R α , with strength λ α . For a system interacting weakly with each bath in the absence of the perturbation δ = 0, assuming the baths are thermal and initially uncorrelated with each other, the reduced state of the system reads is the αth bath correlation function, ρ α is the state of the αth bath, and S α (ω) are the Fourier decomposition of the operator S α which are expressed using the Bohr frequencies of H 0 through where Π are the projection operators onto the energy eigenspace of H 0 . The sums over ω and ω index the system Bohr frequencies. Written in this form, it is clear that the system's eigenstructure is encoded in the propagator. For details see App. A.
Equation (6) is a generic result of second order perturbation theory in the system-bath coupling. To bring it to the LGKS form, two additional approximations are carried out. The first approximation assumes that the integral on the left-hand side samples the function F (τ ) in sufficient accuracy to justify the Fourier transform on the right-hand side, which defines a time independent relaxation rate Γ α (ω) This is the standard Markovian approximation and is valid for long times such that ωt 1. The second assumption is typically a stronger condition than the first, and is referred to as the secular approximation [48]. This approximation is valid when min{t|ω − ω |} 1. The resultant Markovian master equation in the Schrodinger picture reads [49], with the dissipative part and the relations S α (−ω) = S † α (ω) and Γ(−ω) = e −β ω Γ(ω). For Bosonic bath at thermal equilibrium, one can factorize the relaxation rate, Γ(ω) = γ(ω)(n + 1) and Γ(−ω) = γ(ω)n, where n is the Bose-Einstein distribution and γ(ω) is determined by the coupling strength and the density of modes of the bath. Note that the simple additive structure for multiple baths in Eq. (12) results from linearity of the original Hamiltonian and the assumption of no initial correlations between the baths.

III. RESPONSE THEORY WITHIN THE LGKS FRAMEWORK
Now we consider adding a constant perturbation at time t = 0 to the system, which drives the system away from its current NESS to a new steady state. Our treatment of the perturbation δV does not constrain it to be the smallest energy scale in the problem. It is considered small with respect to the system "bare" Hamiltonian H 0 ; however, it is not necessarily small with respect to the system-bath couplings, {λ α }, which is assumed to be weak. An exact LGKS master equation describing the new dynamics would require deriving the master equation that corresponds to the new Hamiltonian H 0 + δV from first principles. In many cases, this global approach is convoluted or even impossible. Here, we present an approach based on a perturbative treatment of the LGKS equation in such a way that is consistent with the approximations in Eqs. (9) and (10), namely, with the secular and Markovian approximations.
The perturbative treatment begins by expanding the eigenvalues and eigenstates of the system Hamiltonian and the perturbation, H 0 + δV , in the basis of H 0 and in orders of δ, Given the structure in Eq. (6), the expansion suggests that the perturbation may influence the master equation in two ways: 1) modifying the structure of the master equation due to changes in the eigenvectors and 2) modifying the Bohr frequencies of the system and the structure of the master equation, due to changes in the eigenvalues. These two modifications can appear either separately or together, depending on the timescales of the problem and the details of the perturbation.

A. Response due to changes in the eigenvectors
To start, we consider how changes to the eigenvectors impact the system's response. We assume that the Bohr frequencies ω do not change as a consequence of the perturbation. This happens anytime the perturbation shifts all the energy levels by a constant, or more commonly, when the perturbation does not affect the eigenvalues at a given order of δ. For example, whenever the eigenvectors of the Hamiltonian H 0 , and the perturbation V are orthogonal, the first order correction of the eigenvalues vanishes, E (1) n = 0. When only the eigenvectors change, the response theory takes on a transitional perturbative form.
Using Eq. (13), the projection operators can be expressed in orders of δ where Π 0 = |ψ n ψ n |, and the expansion immediately follows. By virtue of Eq. (12) and (14), the expansion of the dissipator for each bath in orders of δ, is analogously clear. Here, D j is the δ j order correction to the dynamics, and the zero order term D 0 corresponds to the dissipator of the unperturbed system with the Hamiltonian H 0 . More specifically, the explicit form of the first order corrections is given by and Similar terms in Eq. (16) refers to all terms contributing to first order in δ. The procedure above provides a scheme to treat perturbatively complex Hamiltonians within the LGKS formalism. It is simplified significantly by the constant frequencies entering into the sum, which thus do not alter the underlying Markovian or secular approximations. The ability to write the expansion Eq. (15), admits a quantum response theory for NESS to be developed in an analogous way for a closed quantum system. After the perturbation is turned on, the dynamics of a system in contact with multiple baths follows the master equation, where for each αth bath we can expand the dissipator in a series. The new perturbed state of the system ρ can similarly be expanded in orders of δ as well, where at time t = 0 the NESS is ρ(0) = π 0 . Note that the ρ i corrections are not proper density matrices in themselves and satisfy Tr [ρ j ] = 0 at each order. Inserting Eq. (19) into Eq. (18) and keeping only the first order in δ terms, we havė with the formal solution for the first order correction to the density matrix in δ.
Then, to first order in δ for an arbitrary observable A, where the average · δ is taken with respect to the perturbed state. This implies a response function, where A(τ ) is the observable propagated according to the unperturbed dynamics, is the dynamical map in the Heisenberg picture with the adjoint LGKS generator L † 0 [48]. Note that this result also holds for an arbitrary initial product steady-state, including thermal equilibrium, in which π 0 ∝ exp[−βH 0 ], with the inverse temperature β = k B T and the Boltzmann factor k B . Extensions to second and higher order are straight-forward under the assumption that the eigenvalues do not change.
The response functions derived here are different form those found in the literature, as they involve a term due to the rotation of the dissipator into the perturbed eigenvectors D 1 . The linear response function now generally has two contributions, φ (1) where the first term φ is expected from Kubo theory. The additional term φ (1,2) A originates from noncommutativity of the new Hamiltonian including the perturbation and the system-bath interaction Hamiltonian, and has no classical counterpart.
For closed quantum systems at thermal equilibrium, which follow a unitary evolution, the term φ (1,2) A vanishes and the response function φ (1) A in Eq. (23) reduces to the standard Kubo formula Eq. (1). However, for open quantum systems this is no longer the case. In Sec. III B, we show that this additional contribution can become significant and cannot be neglected.
We point out that a response function which explicitly includes the dissipative part has been introduced before in the literature [28,32,50,51]. However, in these studies the LGKS-based Kubo formula is derived under the assumption that a Dyson expansion of the full evolution exists, without providing a microscopic derivation, or that the perturbation itself is assumed to be represented by a non-Hermitian operator. By contrast, in this study we show that an Hermitian (Hamiltonian) perturbation of an open quantum system, in the limits discussed above, results in an additional contribution φ (1,2) A to the Kubo formula. Moreover, we provide a method to evaluate this term explicitly from first Hamiltonian consideration.

B. Example: Two coupled anharmonic oscillators
In this section, we illustrate how a perturbation that changes only the eigenvectors of the system to first order results in a linear response determined by both the traditional Kubo term and the additional force from the baths. In particular, we investigate the contribution of the different terms of φ has non-negligible contributions. The system we consider consists of two coupled anharmonic oscillators, a and b, subject to a linear external field. Each of the oscillators are coupled to a thermal bath with temperature T a and T b , respectively as shown schematically in Fig. 2a). The anharmonicity is described by a cubic potential δV , which will be considered the perturbation and is treated within the linear response theory.
The unperturbed system Hamiltonian includes the self Hamiltonians of the two harmonic oscillators, a coupling between the modes, and a linear field.
Here, ω denotes the frequency of both oscillators, g their bilinear coupling, and ε the linear potential that both feel. The system-bath interaction Hamiltonians are assumed to be bilinear, with R α = j λ α,j (r α,j + r † α,j ) being a sum of displacement operators of the αth bosonic bath, and λ α,j the characteristic interaction scale.
The joint master equation for the harmonic system including the linear field is derived using a global approach [34], which assumes that the inter-coupling is strong compared to the system-bath couplings. The resultant adjoint propagator is where we define the Bose-Einstein distribution n α = [exp(β α ω ) − 1] −1 for each bosonic bath at inverse temperature β α , and the decay rates γ . The termD α, † 0,+ arises from the linear field and can be derived using the perturbation theory that was introduced above (see App. B for details).
We will consider for concreteness the response of the occupation number of oscillator a to a step perturbation of the form δV , with which adds a cubic potential rendering the system anharmonic expressible as a sum of the new normal modes. Since δ sets the scale of the perturbation we are free to chose Ω ≡ ω − . For this choice, the perturbative treatment implies δ 1. In the absence of the perturbation, the occupation number can be calculated using the master equation, which is valid within an arbitrary steady state. Following the approach described in previous sections, we find that the first order correction to the eigenvalues vanish, and thus the Bohr frequencies remain the same. The eigenvectors however are modified, leading to corrections of the dissipator in the format of Eq. (17) (see App. B for details). To calculate the response function φ (1) a † a of the observable a † a, we first calculate the evolution of the operator a † a(τ ) = e L † 0 τ a † a and then evaluate L † 1 a † a(τ ) . In App. B, we show that the time dependent operator a † a(τ ) can be expressed explicitly as a linear combination of operators at time zero multiplied by time dependent functions.
To study the contribution from the different terms ∝ δ, we use the splitting of the response function in Eq. (24). Direct calculation results in and where both contributions depend on the average Bose distribution,n − = (n a − + n b − )/2. The response functions scale linearly with ε. When no linear field is present, the linear response will vanish completely. Notice that only φ (1,2) a † a depends on the kinetic parameters describing the coupling to the bath, similar to the frenetic contributions arising in classically driven diffusive systems [19,20].
In Fig. 2 we plot the two contributions to the response of the occupation a † a in time, i.e. the integrated response functions for i = 1, 2. In the high temperature limit illustrated in Fig. 2, it is clear that the term involving φ (1,2) a † a has nonnegligible contributions both in the transient and steadystate regimes.
The significance of this term is also manifested in the scaling behavior of the steady-state response of a † a with the average temperatureT = (T a + T b )/2. Integrating the response functions over the interval t = [0, ∞], yields the steady state response. In the high temperature limit, n α − ∝ T α , which implies that the term φ a † a ∝T 2 . This change in the scaling with the average temperature is depicted in Fig. 3, where we plot the total steady state response, δ(Φ a † a ), as a function of the average temperature. AsT increases, the non-Kubo term, Φ (1,2) a † a , dominates the behavior with a different scaling and a different sign of the response emergeing.
We note that at high occupation number of the oscillator, the perturbation theory breaks down as the perturbed state depends explicitly on the level number n. In this case there is an interplay between the magnitude of δ and n. In particular, we wish δ n − 1 2 (see App. B). Interestingly, at zero temperature, when n a − = n b − = 0, the steady-state response is independent of γ − , which characterizes the bath and its coupling to the system.

C. Response due to changes in the eigenvalues
Changes to the eigenvalues E n leads to modifications of the Bohr frequencies, which can be expanded in terms of δ as well, ω = ω (0) + δω (1) + · · · . The procedure introduced for deriving the LGKS master equation can now be repeated using the new perturbed Bohr frequencies and the replacement of the projectors, Π (0) n → Π n to a desired order. Modifications to the Bohr frequencies introduce new timescales that become relevant for the LGKS approximations in Eqs. (9) and (10) to be carried out. The weak system-bath coupling limit [47,52] that leads to the LGKS master equation assumes a time coarse-graining of fast oscillating terms with frequencies ν = ω − ω. This approximation is valid only when the system's relaxation time τ satisfies min{τ |ν|} 1. The perturbation may introduce new frequencies, for example at first order in δ, τ |ν (1) | 1, terms oscillating with a † a ) as a function of the average temperatureT . The parameters are γ− = γ, ω = 100γ, g/ = 90γ, ε/ = 5γ, δ = 0.02 and the temperature difference is ∆T = 50 γ/kB. frequencies ν (1) = δ(ω (1) − ω (1) ) are eliminated. Physically, the coarse-granining is carried out when the width of the spectral line is smaller than the level splitting that is now proportional to δ.
When the Bohr frequencies are modified by the perturbation, expansion of the dissipator in orders of δ does not necessarily exist. The Fourier transform of the bath correlation functions in Eqs. (9) and (12) includes the perturbed Bohr frequencies Γ α (ω (0) ) → Γ α (δω (1) ). This implies that a first-order correction to the eigenvalues, given the coarse-graining discussed above, already results in all orders of the dissipator, making the response of the open system nonlinear. Moreover, in the limit δ → 0 the dissipator will generally be different from its unperturbed form. The reason is the time coarse-graining, which eliminates terms oscillating with frequencies ν (1) that are proportional to δ. This discontinuity is a known phenomena when deriving the LGKS master equation from first principle [53,54]. The observations above suggest when the perturbation is small with respect to the system Hamiltonian H 0 , but large with respect to the system-bath interactions, response theory for open quantum systems cannot be formulated by standard means. The response function can no longer be associated with correlation functions in the unperturbed system, and requires knowledge of the perturbed state of the system. However, when the perturbation is also small with respect to the system-bath coupling a generalized fluctuation-dissipation relation can be derived. When in addition the first order correction to the eigenvectors vanishes, a local master equation, in which the perturbation does not influence the dissipator, can be applied. In this limit, when the perturbation is small with respect to the system-bath coupling, ν (1) τ −1 , following the development in Sec. III A, φ (1) where contributions from a nonlocal dissipator, that pre-viously were expressed through φ (1,2) A enters as a higher order cross term, and thus in this instance should be neglected. So although the response function cannot always be associated with correlation function in the unperturbed system, a perturbative treatment based on Eq. (13) offers a practical method for treating complex perturbations which cannot be analytically solved exactly.

D. Example: Two coupled qubits
In this section, we study an example illustrating how a perturbation that modifies only the eigenvalues to first order, and does not change the eigenvectors, leads to an LGKS master equation that involves higher orders of the perturbation as discussed in the previous section. While the perturbation is assumed to be small compared to the system Hamiltonian, the response of the system is investigated in two different limits. The first assumes the perturbation is small compared to the system-bath couplings, which in terms of timescale can be translated to ν (1) τ −1 , with τ the typical relaxation time of the system. Since the system-bath coupling is already assumed to be weak, the regime of applicability is quite restricted and a local approach can be implemented. The opposing limit, in which the perturbation is assumed to be large with respect to the system-baths coupling, i.e. ν (1) τ −1 , leads to new Bohr frequencies and to the coarse-grained dynamics discussed above.
In this example, the response of the heat current between the thermal baths and a finite size system is studied. In particular, we consider a chain of two qubits A and B that are coupled to two bosonic heat-baths with temperature T a and T b , respectively. A schematic of the system is shown in Fig. 4a). The unperturbed system Hamiltonian, including the qubits and the coupling between them, reads where σ i is a Pauli matrix, ω is the transition frequency of the uncoupled qubits, and g the coupling between them.
The system-bath coupling Hamiltonian is assumed to be bilinear, where R χ = k λ χ,k (r χ,k + r † α,k ) with χ = {a, b} is the weighted displacement operator for the χth bath. In App. C we derive the global master equation which is used to calculate the NESS properties of the system, noting that the eigenvalues of H 0 are ω ± = ω±g/ each with a twofold degeneracy. The resultant adjoint propagator is given by, with where again n χ l = (exp(β χ ω ) − 1) −1 is the Bose-Einstein distribution at inverse temperature β χ for the χth bath. The operators appearing in D χ 0, satisfy the relation [H 0 , χ ± ] = − ω ± χ ± with χ ± being either . Each channel has a decay rate γ .
We consider the response of the heat flow through the system at steady-state, J , to a perturbation of the form with Ω the coupling energy. At first order, this perturbation changes the eigenvalues, leaving the eigenvectors the same. The average heat flow from the a bath in the absence of the perturbation is given by (48) and App. C). At steady state, we find heat flows in two channels with energies ω ± and corresponding rates proportional to γ ± . In the high temperature limit, the heat current J a 0 ∝ ∆T scales as the temperature difference ∆T = T a − T b .
When the perturbation is much smaller than the coupling to the baths, ν (1) τ −1 , a local master equation with system Hamiltonian H 0 + δV in the commutator and a local dissipative part D 0 can be justified. The basis of this approximation is neglecting terms of the order O(δ 2 λ 2 ) and higher in the master equation, and the fact that the first-order correction of the eigenvectors vanishes. Incorporating such terms in the master equation would require introducing higher order corrections in the system-bath coupling for consistency.
The heat current from the a-bath, to first order, is given by J a 0 + δJ a 1 , where J a 0 is given by Eq.(42), and J a 1 = =± D †a 0, V , which is strictly the local contribution of the perturbation. Here we wish to emphasize that locality does not refer to individual qubits, as these are treated globally (see Eq. (42)), but rather to the perturbation itself. At steady-state, Note that J a 1 has the opposite sign has that of J a 0 , and the overall contribution is determined by the sign of δ. At the high temperature limit, J a 0 ∝ ∆T /T , whereas J a 1 ∝ ∆T /T 2 withT = (T a + T b )/2 the average temperature.

FIG. 4.
a) Schematic illustration of the system. b) The steady-state heat currents with and without (black) the perturbation at low temperatures,T = 50 γ/kB, as function of the temperature difference between the baths. The red line corresponds to J a 0 + δJ a 1 , and the blue line to J a of Eq. (46). c) The steady-state heat current response as a function of the average temperatureT for a fixed temperature difference ∆T = 20 γ/kB . In red, the first-order correction δJ a 1 . In blue, the difference, J a −J a 0 , between the new and old steadystate heat current (subtracting Eq. (42) from Eq. (46). The parameters are γ− = γ+ = γ, ω = 10 3 γ, g/ = 2 · 10 2 γ, , and δΩ/ = 50γ.
In the opposite limit where the perturbation is much larger than the coupling to the bath, ν (1) τ −1 , the baths act to thermalize a different system with a different spectrum. To first order in δ the eigenvalues and the corresponding eigenvectors read which in this case also happens to be the exact correction for any δ.
The Fourier decomposition of the interaction Hamiltonian operators can be expressed as , with the corresponding frequencies Following this decomposition a master equation can be derived (see App. C), and the heat current from the a-bath reads An analytic expression of Eq. (46) exists, however, it is to involved to present here. Note that averages in Eq. (46) are taken with respect to the new nonequilibrium state. While for the unperturbed system heat was flowing in two channels, in the current limit, heat flows in four channels that correspond to the perturbed Bohr frequencies {ω j } of Eq. (45). Figure 4b) illustrates the steady-state heat currents with and without the perturbation at low temperatures and as function of the temperature difference between the two baths. The black dashed line corresponds to the heat flow of the unperturbed system J a 0 . The red and blue lines indicate the heat flow subject to perturbation calculated when ν (1) τ −1 and ν (1) τ −1 , respectively. Note that according to Eq. (46) we have ν (1) = 4δΩ/ . When first order corrections to the Bohr frequencies are accounted for, we observe a significant amplification of the steady-state heat current (blue line). Whereas, standard response theory predicts only very small changes to the heat flow.
An informative comparison between the limits is obtained by plotting the response of the heat current at steady-state as a function of the average temperatureT Fig. 4c). The red line corresponds to the first-order correction δJ a 1 , and the blue to the difference between the new and old steady-state heat flow, J a −J a 0 , of Eqs. (42) and (46). At low temperatures, the behavior of the heat current response in the two limits is substantially different. This clearly indicates that the perturbation has a significant effect on the dissipation caused by the coupling to the baths. While in both cases the perturbation is considered small with respect to the system Hamiltonian, its relation to the relaxation rate results in a different response behavior. In the high temperature limit, the correction to the Bohr frequencies becomes negligible, the heat current responses in the two limits coincide and vanish asymptotically, and the relaxation rate increases with the temperature.
The fact that a first order perturbation, with respect to the system Hamiltonian, can lead to a nonlinear response of the heat current is illustrated in Fig. (5). In this figure we plot the steady-state heat current as function of the perturbation strength δ. The black dashed line represent the heat current of the unperturbed system given by Eq. (42), the red line is the sum J a 0 + δJ a 1 of Eqs. (42) and (43) that correspond to the limit ν (1) τ −1 , and the blue line is the heat current given by Eq. (46) in the apposing limit ν (1) τ −1 . In the low temperature limit, small changes in the eigenvalues lead to pronounced and nonlinear changes in the heat current with respect to the unperturbed heat current. In the high temperature limit, as may be expected, the perturbation has very little effect on the heat current, as it is dominated by the temperature.

IV. HEAT FLOWS, ENTROPY AND ENTROPY PRODUCTION
As stressed by Alicki and others [34,36,47,54], although LGKS operators cannot generically be associated with an underling Hamiltonian, a thermodynamic consistency requires this association, as a consistent microscopic description of dissipation is not otherwise possible. To accomplish this faithfully, the LGKS operators requires information on the eigenstrcture of the system whose transitions they promote, as encoded by the detail balance relationship between their transition rates. As we have shown above, to accurately obey these constraints in the presence of a perturbation requires care.
To illustrate the thermodynamic compliance of the theory discussed above, we introduce the response of the heat flows, entropy and entropy production. While Eq. (23) provides the response function of a general quantum observable of the system, in this section we introduce the response of thermodynamic functions which are state dependent. The entropy production σ has contributions from changes in the internal entropy of the system and from the exchange of heat with the baths at different temperatures Here S(t) = −k B Tr [ρ(t) ln ρ(t)] is the von-Neumann entropy and J α is the heat flow from the α-bath The state π α appearing in Eq. (48) is the thermal state of the generator D α , with the temperature β −1 α = k B T α .
The average, on the other hand, is taken with respect to the time dependent state, · ρ = Tr [ρ(t)·]. Noting that using D = α D α , Eqs. (47) and (48) together with Spohn's inequality −Tr [(Lρ)(ln ρ − ln π)] ≥ 0 [55], we conclude that the entropy production is non-negative σ(t) ≥ 0 (see App. D). In case the eigenvalues are perturbed, as discussed in Sec. III C, the modified LGKS generator of each bath has a unique Gibbs-like stationary state that corresponds to the perturbed spectrum. This is demonstrated explicitly in the example of Sec. III D in which the perturbation treatment and the global approach coincide. Following similar arguments discussed above, the perturbed entropy production is nonnegative as well. This is detailed in Apps. C 2 and D 1. We mention in passing that if particles are also exchanged between the baths and the system, one needs to account for the heat carried by them, see for example [56].
Next, we derive the linear response for the case discussed in Sec. III A, when the perturbation changes only the eigenvectors. At time t = 0, before the perturbation is applied, the steady state heat flow from the α-bath is strictly given by J α 0 = −β −1 α D α † 0 ln π α 0 . The first order correction, J α 1 = ∂ δ J α | δ=0 , can be expressed in two useful forms. The first reads where π α ∝ e −βα(H0+δV ) is a unique stationary state of the generator L α , see App. D for details. The second form, shows clearly that the sum of the heat flows from all the baths returns the total energy change in the system. The first order correction of the von-Neumann entropy, allows us to calculate the correction for the change in the internal entropy, see App. D, In the second equality, we applied Eq. (21), and in the last, we integrated using the relations Λ † 0 (τ )L † 0 = ∂ ∂τ Λ † 0 (τ ), and Λ † 0 (0) = 1. Assuming the existence of a new stationary state, one may observe that, as may be expected, for t → ∞ the derivative d dt S 1 = 0. To show this, we assume ρ 1 (∞) = π 1 such that L(π 0 + π 1 ) = O(δ 2 ). Then, L 0 π 1 = −L 1 π 0 , which implies that for sufficiently long times, the first line in Eq. (53) vanishes.
At equilibrium, the first order correction to the entropy production vanishes at all times σ 1 (t) = 0. The change in the von-Neumann entropy equals the entropy flow from the bath. For systems at NESS, this is no longer the case. The correction can turn negative or positive.

V. OUTLOOK
Obtaining an accurate description of open quantum systems' response to perturbation is imperative for understanding the dynamics and manipulating the outcome and performance of quantum devices. In this study, we employed a microscopic Hamiltonian approach based on physical arguments, rather than a phenomenological approach, to develop a quantum response theory. This distinction is crucial, as one cannot associate an arbitrary Hamiltonian with a given dissipator of the dynamics. The perturbation modifies the system and the way it is perceived by the surroundings. As a result, the dissipator is modified, and in the linear response, corrections to the standard response function φ (1,1) A are already expected.
Combining stationary Hamiltonian perturbation theory with the LGKS formalism, we developed a useful practical scheme to obtain the response functions for both systems at equilibrium or NESS. The approach reveals the role of changes in the eigenvaectors and eigenvelues as a result of the perturbation. Eigenvector modifications lead to corrections of the dissipator that can be expressed order by order in the perturbation parameter. It therefore allows us to derive response functions that include new non-negligible terms. For closed quantum systems at thermal equilibrium, our result recovers the standard Kubo formula. The new terms can be interpreted as a response to forces induced by the surroundings as a result of the perturbation.
Changes in the eigenvalues lead to modifications of the Bohr frequencies, and in turn to modifications of the operators appearing in the dissipator. In this scenario, a distinction between two limits is appropriate. First, when the system-bath coupling is small compared to the perturbation, linear changes to the eigenvalues lead to the system's nonlinear response. Although now the response function cannot be expressed by correlation functions at equilibrium or NESS, the perturbative approach introduces a genuine scheme treating complex Hamiltonians within the LGKS formalism.
In the opposing limit, the coarse-graining procedure that leads to the LGKS master equation is not justified. However, when the eigenvectors are not influenced by the perturbation and since the perturbation is a assumed to be small compared to the system-bath coupling, a local master equation is justified and the linear response is strictly given by φ (1,1) A of Eq. (24). Taking into account higher-order contributions would make sense only if higher orders in the system-baths coupling are considered as well. While in this work we focused on stationary perturbations, there is still a need for a quantum response theory of time-dependent perturbations that are derivable from a physical Hamiltonian perspective.
We briefly review the main ingredients for deriving the mater equation in the weak coupling limit. More elaborate derivations can be found in Refs. [48,52,57]. We assume a quantum system with the Hamiltonian H 0 = m E m |m m| coupled to a thermal bath ρ R with the Hamiltonian H R such that [ρ R , H R ] = 0, via an interaction Hamiltonian H I = λS ⊗ R. Here, λ represents the coupling strength which is assumed small, and S and R are linear operators of the system and bath respectively. The generalization to a system coupled to several bath is straight forward and follow the linearization of the master equation for initially uncorrelated baths.
Working in the interaction picture, the reduced dynamics of the system is given by the partial trace over the bath degrees of freedom where is the time ordered propagator in the interaction picture, which is defined using (A3) As shown in [47,53], the dynamical map Λ(t) can be expressed by the cumulant expansion Since we assume the state of the bath is thermal, then Tr [ρ R R] = 0, which implies that the first order cumulant K (1) = 0. The Born approximation (weak coupling) consists of terminating the cumulant expansion at n = 2, henceforth we denote K (2) ≡ K and One obtains with the bath correlations F (s) = Tr(ρ R R(s)R). The Markov approximation (in the interaction picture) means that for long enough time (or short correlation time) one can use the following approximation where L is a LGKS generator. To find its form we first decompose S(t) into its Fourier components, which in the interaction picture reads where the set {ω = E n − E m } contains the Bohr frequencies of the system Hamiltonian. This decomposition is equivalent to the requirement [H 0 , S(ω)] = ωS(ω). Expression (A7) then reads, To bring (A9) to the LGKS form, two approximations are carried out The first approximation assumes that the integral on the left-hand side sample the function F (τ ) in sufficient accuracy in order to justify the Fourier transform on the righthand side. This approximation is valid for long times such that t 1/ω. The second assumption is typically a stronger condition than the first, and is referred to as the secular approximation. This approximation is valid when t max{1/|ω − ω |}. Finally, the Markovian master equation in the Schrodinger picture reads Moving to the interaction picture we have Next, the derivation of the master equation in the weak coupling limit can now be performed as presented in appendix A. When the difference between the new Bohr frequencies |ω + − ω − | = 2g/ is greater than the relaxation rate, terms rotating at that frequency are neglected and the two bosonic modes d + and d − can be considered as independent harmonic oscillators, where each is coupled to both baths. The joint master equation for the Harmonic system including the linear field is then given by, Here we defined the Bose-Einstein distribution n α = [exp( ω β α ) − 1] −1 , and the decay rates γ ≡ γ(ω ), with = ± and α = a, b. The termD α † 0,+ arising from the linear field, can also be derived using the perturbation theory that was introduced above. In this case, a first order perturbation of the linear field recovers the exact master equation, just like a perturbation theory for the harmonic oscillator in a linear field of a closed system.

Cubic perturbation
We study the linear response of the system above to a cubic potential of the form we are left with evaluating Averaging with respect to the NESS we are left with and using the steady state averages we arrive at the result Where In this limit, oscillating terms with frequencies 4δΩ/ vanish and the master equation takes the form with D i j ρ = γ j (n i j + 1) S j ρS † j − + γ j n i j S † j ρS j − 1 2 {S j S † j , ρ} .
In the weak coupling limit, energy and entropy flows from the baths can be associated to changes in the energy of the quantum system determined by the partial generators L α of the α-bath. The change in the energy of a quantum system with Hamiltonian H can be written as Since we assume that L α ≡ −i/ [H, ·]+D α is the thermal generator of the α-bath, each L α as a unique Gibbs like stationary state π α ∝ e −βαH and we can write Putting these together we see that the entropy production is always positive The inequality follows from Spohn inequality −Tr [(Lρ)(ln ρ − ln π)] ≥ 0 that apply to any LGKS generator with a stationary state π. Since any L α has a unique stationary state Π α then also the sum over α in Eq. (D4) is nonnegative.

First order correction of the heat flows
The first order corrections for the heat flow can be derived in two ways. In the first, the starting point is by considering the total energy change of the system. The first order correction reads ∂ δ Tr ρL † H | δ=0 = Tr ρ 1 L † 0 H 0 + Tr π 0 L † 1 H 0 (D5) In the last equality we used eq.(21) for ρ 1 , the relation ∂ τ Λ † 0 H 0 = Λ † 0 L † 0 H 0 , and taking the integral explicitly. One can also note that at steady state, ρ 1 → π 1 , the energy change vanishes as expected. This is an immediate consequence of the relation L 0 π 1 = −L 1 π 0 which holds at steady state. Based on the first equality in (D5) we can identify the correction from the α-bath The second approach is by expanding Eq. (D2), that gives where π α is the new steady state of the generator L α . Assuming the α-bath generator L α has a unique stationary state π α ∝ e −βα(H0+δV ) , equations (D6) and (D7) become identical.