Non-equilibrium quantum probing through linear response

The formalism of linear response theory can be extended to encompass physical situations where an open quantum system evolves towards a non-equilibrium steady-state. Here, we use the framework put forward by Konopik and Lutz [Phys. Rev. Research {\bf 1}, 033156 (2019)] to go beyond unitary perturbations of the dynamics. Considering an open system comprised of two coupled quantum harmonic oscillators, we study the system's response to unitary perturbations, affecting the Hamiltonian dynamics, as well as non-unitary perturbations, affecting the properties of the environment, e.g., its temperature and squeezing. We show that linear response, combined with a quantum probing approach, can effectively provide valuable quantitative information about the perturbation and characteristics of the environment, even in cases of non-unitary dynamics.


I. INTRODUCTION
The study of non-equilibrium scenarios is key for the dynamics -and thermodynamics -of quantum systems [1,2].Although an overarching framework to consistently interpret the variety of non-equilibrium quantum phenomena is hitherto missing, some light has been shed towards understanding several aspects of this multifaceted problem.These theoretical studies come together with a strong interest in applications concerned with transport phenomena, e.g., in quantum technologies at the nanoscale [3].
Given a non-equilibrium setting, the standard approach follows classical statistical mechanics [4], where one applies small perturbations to the system of interest to gather information about the system itself: depending on the nature of the perturbation, one would be able to deduce relevant physical properties through, e.g., response and relaxation functions or generalized susceptibilities [5].This simple, yet insightful, idea is the essence of linear response theory; this being usually the first step in trying to harness the otherwise complex phenomenology of non-equilibrium systems.However, it is worth emphasising that, when extended to quantum systems, Kubo's original formulation of linear response theory relies on two assumptions [6,7]: the system is isolated so that its dynamics are unitary; furthermore, it relaxes towards a thermal equilibrium state.These conditions are not always met in practice.On one hand, closed quantum systems are usually a crude approximation, as we should effectively include environmental effects in the form of dissipation and/or decoherence [8,9].On the other hand, there are cases in which the system relaxes towards a nonequilibrium steady-state (NESS).This class of states is obtained, for instance, when considering boundary-driven systems, i.e., systems that are dissipatively driven by coupling with two different baths at the two ends [10].The physical properties of the system and its bath (e.g., temperature or chemical potential), as well the intra-system and system-bath coupling will ultimately influence the transport properties.
In this work, we aim to probe the non-equilibrium dynamics of a system through linear response theory.To this end, we extend the customary domain of application of linear response theory to include non-unitary perturbations to the system dynamics.Non-unitary perturbations have been considered in Ref. [11], where a method was put forward for finding the linear response of the covariance matrix to a perturbation in the Gaussian channel describing its evolution.The perturbation could be either unitary or non-unitary.Here we proceed via a different method where we find the linear response of an observable of the system to a perturbation in the master equation describing its evolution.As a paradigmatic setting of boundary-driven systems, we consider a simple yet highly non-trivial example and address two coupled harmonic oscillators locally interacting with their own bath.The non-vanishing interaction between the two parties eventually leads the composite system to a NESS.Furthermore, we assume that the dynamics of such a system are described by a local master equation which does not include memory effects, i.e., we rely on a fully Markovian description of the evolution [8].This scenario has been analysed in Ref. [12], where, inspired by considerations coming from classical statistical mechanics [13][14][15], a framework for non-equilibrium quantum linear response has been put forth.However, consistently with the original spirit of the linear response approach, the interest is usually limited to unitary perturbations [12,16] -e.g., to the interaction Hamiltonian between the two subsystems [12].
Here we use one of the two oscillators to test the effects of the perturbations to the bath that is locally coupled to the other one.This configuration is fully in line with the quantum probing paradigm, where it is customary to use a small, simple, and controllable system interacting with a more complex environment to infer accurate information about environmental parameters (e.g., temperature or spectral properties) [17][18][19][20].In turn, the problem that we address is closely connected with quantum estimation theory [21][22][23][24][25], for which the ultimate goal is to find the optimal measurement scheme to gain precise information about a set of given parameters characterising the system and its dynamics [26][27][28][29][30]. Interestingly, within the context of non-equilibrium open quantum systems described by Markovian dynamical maps, a close relationship between fluctuation-dissipation theorems [31] and quantum metrology has been established [32], though limited to the static case: such a connection allows for the description of the effect of a perturbation within the linear response regime.
We first investigate the case where the probe itself is in contact with a thermal bath; despite the thermal noise, we are able to detect the changes to the environment interacting with the main system, when a sudden quench is applied to its temperature.In line with the general aim of quantitatively assessing the effects of a sudden change of the macroscopic parameters characterising the dissipation, we also consider the case where the system interacts with a thermal squeezed bath, while we assume that the probe is perfectly isolated.We show that linear response extended to the case of non-unitary dynamics remains effective in providing key information on the perturbation affecting the system and, when paired with a quantum probing approach, offers a valuable quantitative tool for the characterisation of the features of an environment.
The remainder of the paper is organized as follows.In Sec.II we set out the theory of linear response of open quantum systems to small perturbations.We try to keep our formalism as general as possible, so that both unitary and non-unitary perturbations can be applied to the system dynamics.In Sec.III we describe our first system of interest: a pair of coupled quantum harmonic oscillators, each interacting with local thermal environments.Given the nature of the systems and states that we consider, we can solve the dynamics analytically using the formalism of Gaussian quantum mechanics without the need for numerical approximations.In Sec.III A we apply two simultaneous perturbations -one unitary and one non-unitary.We then focus on the special case of a perturbation of the bath temperature.In III B we show that the second oscillator can serve as a probe of the response to a perturbation on the first oscillator.In Sec.IV, we detach one of the two subsystems (which becomes our probe) from the corresponding local bath, while we let the other subsystem interact with a local squeezed thermal bath.We then study the response of the probe to perturbations on the parameter controlling the squeezing of the bath.In Sec.V, we present our conclusions, along with an outlook for future directions.

II. LINEAR RESPONSE IN OPEN QUANTUM SYSTEMS
In this Section, we will review the formalism of linear response for open quantum systems.In particular, we will explicitly show that the linear response of a generic observable to a perturbation can be written in terms of the steady-state of the system and the observable of interest, written in the Heisenberg picture, where the timeevolution is solely controlled by the unperturbed dynamics [12].Since our aim is to investigate how the system responds whenever we perturb either the unitary or the non-unitary part of the dynamics, we will write the relevant equations in a general form.Thus we will work at the level of superoperators, with the Liouvillian of the dynamics expressed as the sum of a unitary part and a dissipator in the standard Lindbladian form [8,[33][34][35].The dynamics of the system are described by the Markovian master equation where the Liouvillian L can be expanded (to first order) as where L 0 ρ(t) = −i[H 0 , ρ(t)] + Dρ(t) with H 0 the unperturbed Hamiltonian and D the system's dissipator.Also, ϵ = ϵ(t) is the time-dependent parameter controlling the perturbation.Following lines similar to those in Ref. [12], we assume that, for a fixed value of ϵ, the dynamics possess a stationary state ρ ϵ , i.e., Lρ ϵ = 0, that can be expressed as where ρ 0 is the steady-state of the unperturbed dynamics [36], and ρ 1 its linear (in ϵ) correction due to the perturbation.It is worth emphasising that ρ ϵ is not necessarily an equilibrium state.When that is the case and the perturbation is unitary, we recover the standard Kubo formalism of linear response [7].Moreover, in general, the state at some time t can be expressed as where π 1 (t) has to be traceless in order for ρ(t) to be a physical state.By substituting Eq. (4) into Eq.( 1) and integrating over time we arrive at the final result to first order in ϵ [12,37].It follows that the linear response of a generic observable A can be written as where ⟨A⟩ 0 is the unperturbed expectation value and R A (t) = Tr[Ae L0t L 1 ρ 0 ] is the linear response function.
Note that this response function can be expressed as where L * 0 is the dual of the generator L 0 of the unperturbed dynamics, and A H (t) evolves according to the adjoint master equation [38] that encompasses the dual dissipator D * [•] [8].

III. SYSTEM IN CONTACT WITH LOCAL THERMAL BATHS
First, let us consider a system of two coupled quantum harmonic oscillators, each interacting with a local thermal bath, as depicted in Fig. 1.The unperturbed Hamiltonian of this system reads where ω 1 is the frequency of the first harmonic oscillator and ω 2 is that of the second, which is detuned by a quantity δ, i.e., ω 2 ≡ ω 1 +δ.The interaction between the two oscillators, controlled by the coupling strength λ, is modelled so as to preserve the total number of excitations across the closed system.The open dynamics of the unperturbed system are governed by the master equation where each local dissipator reads (j = 1, 2) Each D j [ρ], defined in terms of the local creation and annihilation operators a † j and a j , is the sum of two terms: the first describes the incoherent loss of excitations, while the second describes incoherent pumping.Here, γ is the damping rate (assumed to be the same for both oscillators), while the mean number of excitations for the j-th thermal bath at inverse temperature Given that the Liouvillian of the unperturbed dynamics is quadratic in the creation and annihilation operators, we can employ the methods of Gaussian quantum mechanics to solve the dynamics exactly.This is a valid approach as far as we restrict our considerations to Gaussian states and, as we will do in the following, to perturbations that are also quadratic forms of the operators of the oscillators.Hence we only need to study the evolution of the first and second moments of the dimensionless position and momentum operators . Given the vector of quadratures Y ≡ (x 1 , p 1 , x 2 , p 2 ), and starting from Eq. ( 10), it is straightforward to check that the first moments Ȳi ≡ ⟨Y i ⟩ are damped to zero [39].Therefore, if we assume -without loss of generality -the first moments to be initially zero, the system dynamics will be solely determined by the second moments.This is a customary assumption, e.g., when one considers optomechanical set-ups, where the Gaussian modes always represent zero-mean fluctuations around some semi-classical steady-state [40,41].
The second moments are encoded in the covariance matrix (CM) σ, whose entries are defined as The time-evolution is thus expressed in the Lyapunov form as σ = ασ + σα T + D, where α and D are the drift and diffusion matrices, respectively.While the former bears dependence upon both the unitary and the non-unitary (dissipative) part of the dynamics, the latter only depends on the non-unitary part, describing the effect of thermal noise due to the interaction between the system and the environment.In order to obtain the linear response of an observable A of the system as given in Eq. ( 6), we need the steady-state of Eq. ( 10) and the Heisenberg-evolved observable A H (t). The steady-state ρ 0 is completely characterized by the solution to ασ + σα T = −D, which can be found to read as [cf.Appendix A for the details of the calculations] where with While, in general, the coupling between the two oscillators results in a non-thermal steady-state, in the limit λ → 0, we have B → 0, C → 0, and ζ → 1, thus recovering the case of two independent quantum harmonic oscillators, each relaxing towards the corresponding thermal state.We start by considering the observable A 1 = β 1 ω 1 a † 1 a 1 .The corresponding adjoint master equation [cf.Eq. ( 8)] leads to a set of coupled differential equations, whose solution gives us The explicit forms of the functions f (t), j(t), p(t), q(t), and s(t) are given in Appendix B.
A. Two perturbations: coupling strength and temperature Given the unperturbed dynamics in Eq. ( 10), we investigate the response to two perturbations, L ϵ 1 and L ϕ 1 : the former is a sudden quench in the coupling strength λ between the two harmonic oscillators, while the latter is a quench in the number of excitations n 1 of the first bath.Note that the linear response to a Hamiltonian perturbation of the coupling strength has been thoroughly studied in Ref. [12].The perturbed master equation reads We take L ϵ , where H P is a beamsplitterlike Hamiltonian perturbation of the form [42] where ϵ > 0, and θ(t) is the Heaviside step function.
We obtain the linear response of observable A 1 due to to the Hamiltonian perturbation L ϵ 1 as Since , we can use Eq. ( 14), and rewrite the expectation values in terms of entries of the steady-state CM from Eq. ( 12).We get where we have introduced the kernel function with ∆n = n 2 − n 1 and z 2 = δ 2 + 4λ 2 .Consider now a quench in the temperature of the environment of the first oscillator through a perturbation in n 1 described by where ϕ(t) = ϕ θ(t), with ϕ > 0. The corresponding linear response is After some algebra, we find with Finally, the linear response to a combination of perturbations to both the coupling strength of the interaction and the temperature of the first oscillator's bath is obtained by simply adding the two aforementioned contributions.Therefore, the linear response (of the energy of the first oscillator) to both perturbations is where ∆A ϵ 1 (τ ) and ∆A ϕ 1 (τ ) are given by Eqs. ( 18) and ( 22), respectively.
In Fig. 2, we show this response as a function of time for a certain choice of the parameters of the system (see caption) and with perturbation strengths ϵ = 0.1ω 1 and ϕ = 0.1.We can also compare the linear response result with the exact dynamics of the perturbed system.In order to do so, we consider the steady-state ρ ′ 0 reached asymptotically by the system when setting the coupling strength to λ + ϵ and the average number of excitations in the bath of the first oscillator to n 1 + ϕ, i.e., taking the perturbations into account.The difference between this new steady-state and the original unperturbed steadystate shall be referred to as the perturbed value of the energy response of the oscillator.We first calculate the expectation value of A 1 over the steady-state ρ ′ 0 , i.e., , as we are looking for asymptotic conditions.The invariance of the steadystate under the dynamics allows us to remove the time dependence, thus leaving where (σ ′ ij ) 0 are the elements of the CM associated with ρ ′ 0 .From this, we need to subtract the expectation value calculated over the unperturbed steady-state, for which Using the CM given in Eq. ( 12) we can find the explicit form of the perturbed value of the energy response of the first oscillator, and compare it with the long time limit of Eq. (23).It is straightforward to show that the two expressions coincide to first order in ϵ and ϕ, consistently with the linear response theory approximation.In Fig. 2, we clearly see that the perturbations cause oscillations in the steady-state response.The amplitude of such oscillations is damped over time, until the linear response converges to the perturbed value as t → ∞.We also consider the case of λ → 0 where the oscillators are uncoupled and only interact with their own local bath.As each oscillator is in its own equilibrium state just before they are perturbed, we refer to this as the equilibrium response.Additionally, the case in which the system evolves unitarily can be recovered from Eq. ( 23) by taking γ → 0. As expected, in this limit, we obtain reversible dynamics, so that the linear response displays recurrences over time, which essentially follow the energy flowing back and forth from one oscillator to the other.
We now analyse a special case of the previous scenario, where we perturb only the temperature of the first bath, while we keep the coupling strength constant.Hence, the response is given by Eq. ( 22), or, equivalently, by setting ϵ = 0 in Eq. ( 23).The results are plotted in Fig. 3.In this case, where the perturbation affects only temperature, the linear response fully captures the exact dynamic response, which is calculated in Appendix.D. The equilibrium response here is identical to the equilibrium response in the case of two perturbations, consistently with the fact Eq. ( 18) vanishes for λ → 0, reducing the case of two perturbations to just the response to the quench in temperature.In the equilibrium scenario, we have f (t) = e −γt , thus Eq. ( 22) yields ∆A 0 1 (τ ) = β 1 ω 1 ϕ(1 − e −γτ ), which, as τ → ∞, converges to the asymptotic value ∆A 0 1 (∞) = β 1 ω 1 ϕ.This is the maximum response possible, and it occurs when λ → 0, i.e., when the oscillators are uncoupled.By contrast, whenever we switch on the coupling between the two oscillators, i.e., λ ̸ = 0, the response of the first oscillator is comparatively smaller, as -through the interaction with the second oscillator -the first oscillator additionally experiences the effects of the second bath.

B.
Second oscillator as a probe of a perturbation Let us now consider the effect of a perturbation to the temperature of the first bath on the energy of the second oscillator.In order to do so, we apply a quench to the number of excitations in the first bath, and investigate the possibility of gaining information about this perturbation through the second oscillator, thus probing the environment of one subsystem through the response of the other.The linear response is once more given by Eq. ( 21) with A 2 = β 2 ω 2 a † 2 a 2 replacing A 1 .The response shown in Fig. 4 is qualitatively similar to that in Fig. 3, but the magnitude is smaller due to the fact that β 2 ω 2 < β 1 ω 1 .Furthermore, the perturbation is now mediated by the first harmonic oscillator, which is the one whose bath is perturbed.As can be deduced from Fig. 4, whenever the coupling between the two oscillators is nonzero, the perturbation applied to the first oscillator affects the second one.As a result, the value of the coupling strength λ determines the response of the second oscillator.We provide evidence of this by exploring the limit λ → 0, which corresponds to the equilibrium response.As we discussed in Sec.III, the equilibrium response of the first oscillator, shown in Fig. 3, approaches the maximum value ∆A 0 1 (∞) = β 1 ω 1 ϕ, whilst the equilibrium response of second oscillator is zero [cf.Fig. 4].For a finite value of the coupling constant λ, i.e., when we look at the steady-state response, the perturbation non-trivially affects the system response.
The difference between the perturbed value for coupled and uncoupled oscillators in Fig. 3 is related to the perturbed value of the second oscillator in Fig. 4 as where ∆A λ j (∞) is the perturbed value of the energy of the j th oscillator (j = 1, 2) for a coupling strength λ.As ∆A 0 1 (∞) = β 1 ω 1 ϕ, we find This expression allows us to calculate the perturbed value of one of the oscillators if we know that of the other.This means that we can use the second oscillator as a probe to gain information about the response to a perturbation to the bath of the first.The degree to which the second oscillator mimics the first is determined by the coupling strength λ.We can add an additional line to the plots in Fig. 3 and Fig. 4 to show the behaviour as λ − → ∞ , illustrated in Fig. 5. Comparing Fig. 5 (a) and (b), we see that when λ − → ∞ , the energy response of the first oscillator approaches a value of β 1 ω 1 ϕ/2, while -in the same limit -the response of the second oscillator tends towards β 2 ω 2 ϕ/2.This shows that in the limit of infinite coupling strength, the energy of both oscillators increases by a factor proportional to ϕ/2.This can be interpreted as an instance of equipartition of energy between the two parties of the composite system.

IV. SYSTEM INTERACTING WITH A SINGLE SQUEEZED THERMAL BATH
We now investigate a system comprised of two coupled quantum harmonic oscillators, one of which is connected to a squeezed thermal bath, as shown in Fig. 6.The environment is characterized by two parameters, i.e., its temperature T and the squeezing parameter s, which can be written in polar form as s = re iθ .The unitary dynamics of this system are described by the Hamiltonian H 0 given by Eq. ( 9).In this scenario, the unperturbed dynamics are governed by the following master equation [43] , 6. Sketch of the second scenario: two coupled quantum harmonic oscillators, one of which is interacting with a squeezed thermal bath characterized by temperature T and squeezing s.We investigate the response of the system directly interacting with the bath through the second system, which serves as a probe, when applying a sudden quench to parameter M , i.e., M → M + η.
where the dissipator -acting only on the first systemcomprises the superoperators with N = n cosh 2 r + sinh 2 r + sinh 2 r, M = − cosh r sinh re iθ (2n + 1), and n = [exp(βω 1 ) − 1] −1 is the average number excitations in the bath at inverse temperature β = T −1 .The parameters M and N are not mutually independent as the condition |M | 2 ≤ N (N + 1) must be enforced in order to ensure positivity of the density matrix [8].
As in the case discussed in Sec.III, we need to find the steady-state ρ 0 or, equivalently, the corresponding CM.To do so, we follow a slightly different approach compared to that discussed in Sec.III.The core idea is to remapresorting to a standard set of correspondence rules [44] the master equation from Eq. ( 27) into a Fokker-Planck equation for the characteristic function of our two-mode system, defined as where is the displacement operator associated to the i-th mode [43,45].The steadystate characteristic function χ 0 (α 1 , α 2 ) is then found by imposing the condition χ = 0.By inspection, from χ 0 (α 1 , α 2 ) one can construct the steady-state CM, i.e., Σ 0 , written in terms of the annihilation and creation op-erators.In general, the entries of Σ are defined as with X = (a 1 , a † 1 , a 2 , a † 2 ).Note that Σ and σ are related by the transformation σ = ΛΣΛ † , where Λ is defined as Λ ≡ i=1,2 Λ i with The details of the derivation of the steady-state CM are given in the Appendix C, whereas here we state the final result, i.e., where the explicit form of the coefficients D, E, and F is provided in Appendix C.These coefficients ultimately depend on the physical parameters characterising the dynamics of the system.In particular, in the limit λ → 0, both E and F vanish, as one can immediately check through Eqs.(C6) and (C7).In other terms, when the two oscillators are uncoupled, the steady-state is given by a block-diagonal CM in the form where the steady-state CM associated with the first mode, i.e., is non-diagonal due to the squeezing of the associated bath, whereas the second mode converges to the diagonal CM Σ D = diag (N + 1/2, N + 1/2).Furthermore, in the limit M → 0, we have that the coefficient D vanishes and the system relaxes towards a global thermal state.
A. Perturbation to bath squeezing Given this system, we apply the linear response formalism outlined in Sec.II to study the effect of a perturbation to the squeezing of the first bath.More concretely, the perturbation is described by the superoperator We restrict our attention to the case of a sudden perturbation, which we model through the step-like function The link between M and N entails that a step-like perturbation on the former will in turn cause a similar perturbation on the latter of the same form as L ϕ 1 from Eq. ( 20).The values of η and ϕ are related as under the assumption that the squeezing parameter s = re iθ is not directly affected by the perturbation.The linear response is obtained by accounting for both the perturbations in Eq. ( 20) and Eq. ( 36), i.e.
As we aim to use the second oscillator as a probe, we choose the observable A 2 (t) = (a † 2 a 2 )(t) to determine the linear response in dimensionless units.The calculation of the time-evolution of A 2 (t) for this situation follows lines analogous to those previously illustrated, leading to where f (t), g(t), j(t), and l(t) are defined in Appendix B. Plugging Eq. (39) into Eq.( 38), we find that the contribution coming from L η 1 is zero, hence the linear response reduces to which is of the same form as that used in Sec.III, provided that the time-evolution of A 2 (τ − t) is obtained from Eq. ( 39), and the steady-state ρ 0 is deduced from the CM in Eq. (33).The linear response is finally given by which is plotted in Fig. 7.
The perturbed value of the energy of the second oscillator is given by ϕ = ∆N , where the change in N is caused by the perturbation of the parameter M .At long times, the measurement of the change in the energy of the second oscillator gives a value for ϕ.Hence we can deduce the change in squeezing, η = ∆M , of the bath attached to the first oscillator from Eq. (37).This shows that we can use the second oscillator as a probe of the perturbation to the squeezing of the bath interacting with the first oscillator.

V. CONCLUSIONS AND OUTLOOK
We have used the formalism of linear response theory to investigate both unitary and non-unitary perturbations around non-equilibrium steady-states of an open quantum system consisting of two coupled harmonic oscillators.We have calculated the response to simultaneous perturbations to the coupling strength between the two oscillators (i.e., a unitary perturbation), and to the temperature of the thermal bath (i.e., a non-unitary perturbation).We also looked at the effect solely determined by the temperature perturbation, investigating the response of both the oscillators to this.We found that such an approach is effective in probing perturbations to the dynamics of the system, including in situations where the latter is affected with a non-equilibrium bath.
The approach presented in this paper shows how to extend the theory of quantum linear response to nonunitary perturbations affecting the dynamics of open systems.Our study calls for a mathematically consistent formalism able to microscopically account for such perturbations.This would allow quantitative analysis of more general scenarios where the dynamical equations go beyond the usual Markovian approximation, encompassing, e.g., strong-coupling and memory effects.
Appendix A: Steady-state covariance matrix for the first case For the system described in Sec.III, we find the steadystate using the method set out in Appendix C of Ref. [11] to transform some classes of Lindbladian master equations into dynamical equations for the first and second moments.Let us consider a Lindbladian master equation in the form where the Hamiltonian H 0 is quadratic in the quadrature operators, i.e., it can be written in the following matrix form where Y is the vector of quadratures, while the Lindbladian operators are linear, i.e., L k = c T k Y T .It is immediate to see that the master equation given by Eq. ( 10) fulfils these requirements, provided that we rewrite the relevant quantities in terms of the quadratures x i and p i , which are related to the set of creation and annihilation operators, a † i and a i , by a linear transformation, as outlined in Sec.III.Therefore, the free Hamiltonian of Eq. ( 9) reads The latter can be brought in the form of Eq. (A2) by taking Y ≡ (x 1 , p 1 , x 2 , p 2 ) as vector of the quadratures, and Comparing the unperturbed master equation in Eq. ( 11) with the general form in terms of the jump operators, i.e., Eq. (A1), we can identify the jump operators as Following Ref. [11], we seek the vectors c k such that L k = c T k Y T ; they are given by The vectors c k allow us to construct the matrix The drift matrix α is obtained via [11] α where Im(M) denotes the imaginary part of a matrix M, and Ω is the symplectic matrix given by Ω = Ω 1 ⊕ Ω 2 , where Substituting Eqs.(A4) and (A7) into Eq.(A8) yields The matrix D is given by [11] where Re(M) denotes the real part of a matrix M. Substituting Eq. (A7) into (A11), the matrix D turns out to be Once the matrices α and D have been obtained, we can determine the steady-state CM, σ 0 , by solving ασ + σα T + D = 0, which gives where , and . This matrix can be easily brought in the more compact form of Eq. ( 12).

(C7)
Appendix D: Exact dynamic response In the case of applying one perturbation to temperature or squeezing, the linear response, ∆A(τ ), is identical to the exact dynamic response, which we shall denote ∆A E (τ ).We calculate the exact dynamic response in a similar way to the calculation for the perturbed value [cf.Eq. (24], but this time looking for the change in energy as a function of time rather than at the steady-state.For instance, to find the exact dynamic response of the energy of the first oscillator to a perturbation in temperature, we choose A = β 1 ω 1 a † 1 a 1 and calculate In the expansion of a † 1 a 1 (t) given in Eq. (B6), s(t) is the only term that depends on n 1 .Hence, Eq. (D1) simplifies to ∆A E (τ ) = β 1 ω 1 [s(τ, n 1 + ϕ) − s(τ, n 1 )] , (D2) which exactly matches the linear response plotted in Fig. 3.
Similarly, for the perturbation in squeezing, note that l(t) is the only term that depends on the perturbed parameter, i.e., n, in the expansion of our observed quantity [cf.Eq. (B8)], which is A = a † 2 a 2 in this case.Therefore, the exact dynamic response of the second oscillator to a perturbation in the squeezing of the first bath is given by ∆A E (τ ) = l(τ, n+ϕ)− l(τ, n).Again, the linear response perfectly reproduces this exact dynamic response.

FIG. 1 .
FIG.1.Sketch of the system investigated in Sec.III.It consists of two coupled quantum harmonic oscillators, interacting with their local thermal baths at two different temperatures, T1 and T2.We study the response of one of the two subsystems when we apply a sudden quench either to the coupling constant (i.e., λ → λ + ϵ) or to the number of excitations in the first bath (i.e., n1 → n1 + ϕ).

FIG. 2 .
FIG. 2. Plot of the dimensionless energy response of the first oscillator (with A1 = β1ω1a † 1 a1) to a step perturbation to both the coupling strength λ, and the number n1 of excitations of the first bath.The dynamics are simulated for the following values of the parameters: δ = 10ω1, λ = 5ω1, γ = 0.5ω1, ϵ = 0.1ω1, ϕ = 0.1, β1ω1 = 0.1, and β2ω1 = 0.001.The steady-state response converges to the perturbed value.The equilibrium response shows what happens in the limit of decoupled oscillators (λ → 0), while the unitary response (γ → 0) reflects the case where the system is closed.

FIG. 3 .
FIG.3.Plot of the linear response of the dimensionless energy of the first oscillator (with A1 = β1ω1a † 1 a1), for the case of a step perturbation to the number of excitations n1 of the bath locally interacting with it.The plots are obtained with the following choice of the various parameters: δ = 10ω1, λ = 5ω1, γ = 0.5ω1, ϕ = 0.1, β1ω1 = 0.1, and β2ω1 = 0.001.The steady-state response approaches the coupled perturbed value (λ ̸ = 0), whilst the equilibrium response (λ → 0) approaches its own perturbed value for the case where the oscillators are not coupled.

FIG. 7 .
FIG. 7. Plot of the dimensionless energy response of the second oscillator to a step perturbation in the squeezing (M ) of the bath attached to the first oscillator.Curves are obtained for the following choice of the physical parameters: δ = 10ω1, λ = 5ω1, γ = 0.5ω1, ϕ = 0.1, βω1 = 0.1.The steady-state response converges to the perturbed value.The equilibrium response shows that the energy of the second oscillator is unchanged when the oscillators are not coupled (i.e., when λ → 0).