Fundamental limits for cooling of linear quantum refrigerators

# Fundamental limits for cooling of linear quantum refrigerators

Nahuel Freitas Departamento de Física, FCEyN, UBA, Pabellón 1, Ciudad Universitaria, 1428 Buenos Aires, Argentina Instituto de Física de Buenos Aires, UBA CONICET, Pabellón 1, Ciudad Universitaria, 1428 Buenos Aires, Argentina    Juan Pablo Paz Departamento de Física, FCEyN, UBA, Pabellón 1, Ciudad Universitaria, 1428 Buenos Aires, Argentina Instituto de Física de Buenos Aires, UBA CONICET, Pabellón 1, Ciudad Universitaria, 1428 Buenos Aires, Argentina
September 21, 2019
###### Abstract

We study the asymptotic dynamics of arbitrary linear quantum open systems which are periodically driven while coupled with generic bosonic reservoirs. We obtain exact results for the heat flowing into the network, which are valid beyond the usual weak coupling or Markovian approximations. We prove the validity of the dynamical third law of thermodynamics (Nernst unattainability principle), showing that the ultimate limit for cooling is imposed by a fundamental heating mechanism which becomes dominant at low temperatures: the non resonant creation of pairs of excitations in the reservoirs induced by the driving field. This quantum effect, which is missed in the weak coupling approximation, restores the unattainability principle whose validity was recently challenged.

Introduction: Is there a fundamental limit for cooling? The third law of thermodynamics formulated by Nernst nernst1926 (the unattainability principle) establishes that it is not possible to reach zero temperature in finite time. But this is not a sacred principle and, in fact, its validity has been challenged kolavr2012 ; levy2012 ; cleuren2012 on the basis of new results obtained in the emergent field of “quantum thermodynamics” kosloff2013 ; anders2015 ; brandao2015 . Challenging macroscopic thermodynamical laws is not a luxury but a necessity when studying energy exchange in a regime dominated by quantum fluctuations. Using modern quantum technologies, systems evolving in these conditions (often out of equilibrium and far from the thermodynamical limit), are now regularly built and manipulated abah2012 ; an2015 ; pekola2015 . Thus, these technologies enabled a new generation of machines, operating as engines or refrigerators in the quantum domain. In fact, to understand devices such as the single ion engineabah2012 , the quantum absorption refrigeratorslevy2012refri , the phononic coolersarrachea2012 , or the quantum thermal transistorjoulain2016 , quantum thermodynamics is essential.

In this Letter we show that, contrary to previous claims kolavr2012 ; cleuren2012 , the unattainability principle is indeed valid for a wide class of quantum refrigerators which includes the most prominent cases mentioned above (phononic and ionic coolers, for example). Most importantly, we show that the ultimate limit for cooling is imposed by a fundamental quantum process, which is next-to-leading order in the coupling between the system and the reservoirs and becomes dominant at sufficiently low temperatures. It consists of the non resonant creation of excitation pairs in the reservoirs, which is closely related with the dynamical Casimir effect (DCE)wilson2011 ; benenti2015 . As this process is not of leading order in the coupling between the system and the reservoirs, its fundamental role was missed in all previous works that were based on a weak coupling approximation. On the basis of those analysis, potential violations of the third law were suggested kolavr2012 ; cleuren2012 . Instead, the unavoidable heating process we characterize enforces the third law.

We present exact results for quantum refrigerators formed by networks of oscillators whose frequencies and couplings are driven by external sources. Their Hamiltonian is ( and are vectors whose components are position and momentum operators satisfying canonical commutation relations, being and positive and symmetric matrices. The network is coupled with arbitrary bosonic reservoirs and the total Hamiltonian is . Each reservoir is formed by independent oscillators whose positions (momenta) are (). Thus, and the interaction is described by the bilinear Hamiltonian (with coupling constants ).

Work and heat, the two central thermodynamical concepts are naturally defined once we notice that the expectation value of the energy of , , changes for two reasons. The variation induced by the external driving is associated with the work , while the change induced by the coupling with each reservoir is associated with heat . Thus, we can write , where work and heat rates are, respectively, and (both work and heat are positive when energy flows into the network).

We restrict to consider a periodic driving for which there is a stable stationary regime (see below). In this regime, the quantum state of is periodic and the average of vanishes. Then, the asymptotic value of the average work in a cycle, defined as and the analogously defined averaged heat rate , satisfy the identity . This is nothing but the first law of thermodynamics for a cyclic process. Notice that heat could be defined differently: Thus, the rate of change of the energy of is . As the interaction does not commute with the total Hamiltonian , the energy lost by is not equal to the one absorbed by (and therefore ), However, as shown in SM B, in the stationary regime, in average, no energy is stored in the interaction terms and the identity holds.

The solution: We sketch the method we use and leave most details for the Supplementary Material (SM). The effect of the environments on the network depends on the so called spectral density where, for each reservoir, we define the matrix . For example, the effect of the environment on the retarded Green function , that determines the response of the network to an impulse and is a central ingredient in our analysis, shows up in its evolution equation which reads as , where a dot denotes the derivative with respect to the first temporal argument. The reservoirs determine the damping kernel which is and the renormalized potential .

We will make the following assumptions: a) The external driving is periodic with period (i.e. ) and that is such that the driving does not induce instabilities via parametric resonance. b) The coupling with the reservoirs is such that the Green function decays exponentially with . In such conditions there is a stationary regime, where all relevant quantities are also -periodic, and Floquet theory enables the exact solution of the problem. Details are shown in the Supplementary Material (SM) A, but the crucial observation is that, for a periodic driving, the function is also -periodic for large and therefore can be expanded as . The Fourier coefficients satisfy a set of equations obtained directly from the one for , which read:

 ^g(i(ω+kωd))−1Ak(ω)+∑j≠kVjAk−j(ω)=1δk,0. (1)

Here, , is the Laplace transform of the Green function of the undriven network and is the transform of (see arrachea2012 ; arrachea2013 for a related approach). Finally, we assume: c) that the initial state has no correlations between and and that each reservoir is initially in a thermal state with temperature . Then, independently of the initial conditions, the asymptotic state of has a Gaussian density matrix that is fully characterized by a -periodic covariance matrix whose coefficients admit a Fourier series expansion. For example, position-momentum correlation function, defined as can be expressed as . Here, the Fourier components can be written, as shown in SM A, as:

 σxpj,k=12∫∞0dω(ω+kωd)Aj(ω)~ν(ω)A†k(ω), (2)

where is the Fourier transform of the so called noise kernel (.

Heat rate: In the asymptotic regime, is the temporal average of , where is a projector over the sites of the network in contact with (see SM B). Using this, we find an important exact result: is the sum of three terms

 ˙¯Qα=˙¯QRPα+˙¯QRHα+˙¯QNRHα. (3)

The first term, , is responsible for the resonant pumping (RP) of energy from (or into) . It reads:

 ˙¯QRPα=∑β≠α∑k∫∞0′dωωp(k)β,α(ω)Nα(ω)−−∑β≠α∑k∫∞0′dω(ω+kωd)pkα,β(ω)Nβ(ω) (4)

where is the Planck distribution and is a positive number, proportional to the probability for the network to couple the mode with frequency in with the one with frequency in . The first term in Eq. (4) is positive and accounts for energy flowing out of : a quantum of energy is lost in and excites the mode in after absorbing energy from the driving. The second term corresponds to the opposite effect: A quantum of energy is lost from and dumped into the mode in after absorbing energy from the driving.

The second term in Eq. (3) is responsible for the resonant heating (RH) of and reads:

 ˙¯QRHα=−∑β≠α∑k∫∞0′dωkωdpkα,α(ω)Nα(ω). (5)

Its interpretation is analogous to the previous one except for the fact that in this case energy is transferred, through the network, between modes and in . In all cases can gain or loose energy , depending on the sign of . However, when monotonically decreases with (as it does for the Planck distribution) the upwards flow of energy is more probable than the downwards one and heats up. A subtlety should be noticed: The lower limit of the frequency integral in both resonant terms (RP and RH) is not but . This is because processes with negative , which correspond to the emission into the driving, can only take place if . The role of the low frequency modes (with ) is crucial, as we now discuss.

The term corresponds to a non resonant heating (NRH) effect and reads:

 ˙¯QNRHα=−∑k>0∫kωd0dωkωdp−kα,α(ω)(Nα(ω)+12)−−∑β≠α∑k>0∫kωd0dω(kωd−ω)p−kα,β(ω)(Nβ(ω)+12)−−∑β≠α∑k>0∫kωd0dωωp−kβ,α(ω)(Nα(ω)+12) (6)

(strictly, this formula is valid if the driving is time reversal invariant, i.e, if ; the general expression is included in SM).

The physical process that gives rise to is rather different from the resonant ones discussed above. In the non resonant case a pair of excitations is created. One of them has energy while the other one has energy (these values add up to the driving energy , notice that only enters in the above expression). As opposed to the resonant case, in NRH, excitations are not transferred between modes but created in pairs from the driving. The three terms in Eq, (6) respectively correspond to the following three cases: i) both excitations are created in ; ii) mode is excited in and mode in ; iii) mode is excited in while mode is excited in . In all cases gains energy: in the first case the net gain is while in the last two gains, respectively, and . As , heats up. Noticeably, the non resonant heating term is the only one surviving when all reservoirs are at zero temperature (in that case both and vanish).

In the absence of driving we have and only survives. In this case, only contributes to Eq. (4) and the transition probabilities are symmetric (i.e., ). Then, when is the hottest reservoir (accordingly, for the coldest one). This is nothing but Kelvin’s version of the second law (i.e., “heat flows from hot to cold” martinez2013 ). In the driven case, Eq (3) can be used to derive the most general form of the second law, that reads . We omit the proof here (but include it in SM C for completeness) and focus on the derivation of the third law.

Cooling and the third law: Since both and , the only way to pump energy out of (and therefore cool it) is to design a process such that . In fact, to pump more heat out of than the one flowing into it, requires us to impose spatial (or temporal) asymmetries to the driving or to the coupling with the reservoirs levy2012 ; ticozzi2014 ; arrachea2012 . However, to prove the validity of the third law it is not necessary to go into the details of cooling processes. For this, one should simply notice that both resonant terms vanish at zero temperature while the non resonant heating still survives. This alone is enough to establish the validity of Nernst unattainability principle: heating dominates at sufficiently low temperatures. Moreover, this also implies that for any cooling process there is a minimum achievable temperature which can be estimated as the one for which the cooling term in becomes comparable to the . The minimal temperature is not universal and its properties for various cooling mechanisms will be analyzed elsewhere. Here we focus on analyzing an important example. We estimate the minimal temperature that can achieved by the cooling strategy proposed in kolavr2012 to violate the Nernst unattainability principle.

To gain some intuition on the behavior of the heat rates we use the following simplifying assumptions: a) We assume that the driving is weak (i.e. if ). Then, using perturbation theory we find that , for . b) We assume that the coupling with the environments is also weak. Then, (the Green function of the undriven network) can be expressed as a sum over the normal modes of the isolated network, whose eigenfrequencies we denote as . Using this, the frequency integrals in Eqs. (4) and (5) can be performed since the integrand is strongly peaked around the eigenfrequencies and their sidebands (the peaks in arise because of their dependence on ). In this way, as shown in SM D) the resonant heat rates can be expressed as a sum over all normal modes and sidebands. Finally, in the low temperature limit, this sum is dominated by the term with the lowest frequency () because the Planck distribution enforces a natural cutoff. c) To simplify the analysis we will consider all reservoirs at the same temperature (the most favorable condition for cooling), and use the harmonic driving (for which only appear in the weak driving limit). In this case we obtain

 ˙¯QRPα=e−Ω0−ωdT0|V01|2(π2/8)Γ0Ω20(Ω20−(Ω0−ωd)2)2∑β≠αI0β(Ω0)I0α(Ω0−ωd)×⎧⎨⎩(Ω0−ωd)−Ω0I0α(Ω0)I0β(Ω0−ωd)I0β(Ω0)I0α(Ω0−ωd)⎫⎬⎭ (7)

where is the decay width of the mode (see below) and denotes the matrix element of in the normal mode with frequency . The above equation is revealing: when the spectral densities are identical the heat rate is negative and absorbs energy. However, if the condition is satisfied and the reservoir loses energy. This cooling condition simply states that the cooling process (i.e., extracting energy from and dumping it in the mode of ) has a higher rate than the heating process (taking energy from and dumpling it in mode in ). The reduction in the heating rate arises because the density of final states is small. As explained in SM D, the same codition implies that .

Even if the cooling condition is satisfied, the heat rate rapidly decreases with the temperature. This is because the thermal factor appearing Eq. (7) decreases with temperature faster than any power law. However, there is an interesting strategy that we could use to maximize the heat rate. For this, as suggested in kolavr2012 , we can use an adaptive method by slowly adjusting the driving frequency in such a way that . In this case, for , the heat rate is

 ˙¯QRPα=π28e|V01|2Ω60T0I0α(T0)∑β≠αI0β(Ω0)Γ0 (8)

It is clear that the ratio is of zero–th order in the coupling strength between the system and the environment. Therefore, the above identity shows that is of first order in the coupling strength. In contrast, it can be seen from Eq. (6) that is of second order in the coupling strength for (since in that case the integration domain does not include any resonance peak). We will use these results to estimate the minimal temperature that can be achieved by this cooling protocol. For this, we must study cooling as a dynamical process.

Only a finite reservoir can be cooled. Thus, we assume that the has a finite heat capacity . As discussed in kolavr2012 , depends on the dimensionality () of the reservoir and scales with temperature as . If the rate at which energy flows away from is sufficiently small, one can think that the environment has a time dependent temperature that satisfies the equation . To solve this equation we need an expression for the heat rate. We could use Eq. (8) provided that the rate of change of is smaller than the driving frequency (since in that case we can instantaneously satisfy the adaptive condition ). In this way, we obtain that satisfies

 dTαdt=−1Cv˙Qα∝−γ0η T1+λα−dα (9)

where we used that for low frequencies the spectral density (where is a relaxation rate scaling quadratically with the coupling constants between the system and the environment while the exponent characterizes the environment, being the one corresponding to a ohmic reservoir). Above, is a constant depending on and . The solutions of this equation approach in finite time when , which leads us to the surprising conclusion that the unattainability principle could be violated. However, this argument, presented in kolavr2012 , is not correct, because the above equation for stops being valid at suffiently low temperatures. In that case, , which was not taken into account so far, becomes dominant.

The non resonant heating cannot be neglected in spite of the fact that (for ) it is proportional to . This scaling with makes this term invisible to any treatment based on the weak coupling limit (such as the master equation used in kolavr2012 , which is valid to first order in ). On the contrary, the weak coupling limit captures the resonant term given in Eq. (8), which is first order in (in fact, such expression is equivalent to the heat rate used in kolavr2012 ).

Our analysis, which is non perturbative, shows that the non resonant term given in Eq. (6) will end up stopping any cooling. Moreover, it enables us to estimate the minimum achievable temperature by estimating when resonant and non resonant contributions become comparable. (and is roughly independent of for sufficiently low temperatures) and Eq. (8) shows that the cooling term scales as . Therefore, both terms become comparable for temperatures scaling as . In Figure 1 we see that this naive scaling argument is confirmed by a detailed numerical evaluation of both resonant and non resonant heat rates (the minimal temperature is estimated for various reservoirs characterized by different values of and ploted as a function of ). Thus, the dynamical third law (Nernst unattainability principle) has been restored.

Summary: We presented a complete description of linear quantum refrigerators and demonstrated the validity of the third law. Assuming linearity (i.e. a quadratic Hamiltonian) is not esential for our analysis about the validity of such law. In fact, for very low temperatures the discrete lower sector of the spectrum of the system is the only relevant piece of information. Therefore, our proof of the third law is general. The physical process enforcing the third law is the excitation of the enviromental ground state by the driving, that ends up creating excitation pairs in the reservoirs, as in the dynamical Casimir effect (whose potential relevance for heating, in a different setting, was recently pointed out in benenti2015 ).

This work is supported by grants from Ubacyt, Anpcyt and Conicet (Argentina).

## References

• [1] Walther Nernst. The new heat theorem. Methuen and Company, Ltd, 1926.
• [2] Michal Kolář, David Gelbwaser-Klimovsky, Robert Alicki, and Gershon Kurizki. Quantum bath refrigeration towards absolute zero: Challenging the unattainability principle. Physical review letters, 109(9):090601, 2012.
• [3] Amikam Levy, Robert Alicki, and Ronnie Kosloff. Quantum refrigerators and the third law of thermodynamics. Physical Review E, 85(6):061126, 2012.
• [4] Bart Cleuren, Bob Rutten, and Christian Van den Broeck. Cooling by heating: Refrigeration powered by photons. Physical review letters, 108(12):120603, 2012.
• [5] Ronnie Kosloff. Quantum thermodynamics: a dynamical viewpoint. Entropy, 15(6):2100–2128, 2013.
• [6] Sai Vinjanampathy and Janet Anders. Quantum thermodynamics. arXiv preprint arXiv:1508.06099, 2015.
• [7] Fernando Brandão, Michał Horodecki, Nelly Ng, Jonathan Oppenheim, and Stephanie Wehner. The second laws of quantum thermodynamics. Proceedings of the National Academy of Sciences, 112(11):3275–3279, 2015.
• [8] Obinna Abah, Johannes Rossnagel, Georg Jacob, Sebastian Deffner, Ferdinand Schmidt-Kaler, Kilian Singer, and Eric Lutz. Single-ion heat engine at maximum power. Physical review letters, 109(20):203006, 2012.
• [9] Shuoming An, Jing-Ning Zhang, Mark Um, Dingshun Lv, Yao Lu, Junhua Zhang, Zhang-Qi Yin, HT Quan, and Kihwan Kim. Experimental test of the quantum jarzynski equality with a trapped-ion system. Nature Physics, 11(2):193–199, 2015.
• [10] Jukka P Pekola. Towards quantum thermodynamics in electronic circuits. Nature Physics, 11(2):118–123, 2015.
• [11] Amikam Levy and Ronnie Kosloff. Quantum absorption refrigerator. Physical review letters, 108(7):070604, 2012.
• [12] Liliana Arrachea, Eduardo R Mucciolo, Claudio Chamon, and Rodrigo B Capaz. Microscopic model of a phononic refrigerator. Physical Review B, 86(12):125424, 2012.
• [13] Karl Joulain, Jérémie Drevillon, Younès Ezzahri, and Jose Ordonez-Miranda. Quantum thermal transistor. Physical review letters, 116(20):200601, 2016.
• [14] CM Wilson, Göran Johansson, Arsalan Pourkabirian, Michael Simoen, JR Johansson, Tim Duty, F Nori, and Per Delsing. Observation of the dynamical casimir effect in a superconducting circuit. Nature, 479(7373):376–379, 2011.
• [15] Giuliano Benenti and Giuliano Strini. Dynamical casimir effect and minimal temperature in quantum thermodynamics. Physical Review A, 91(2):020502, 2015.
• [16] Liliana Arrachea and Bruno Rizzo. Nonequilibrium green’s functions in the study of heat transport of driven nanomechanical systems. In Journal of Physics: Conference Series, volume 427, page 012012. IOP Publishing, 2013.
• [17] Esteban A Martinez and Juan Pablo Paz. Dynamics and thermodynamics of linear quantum open systems. Physical review letters, 110(13):130406, 2013.
• [18] Francesco Ticozzi and Lorenza Viola. Quantum resources for purification and cooling: fundamental limits and opportunities. Scientific reports, 4, 2014.
• [19] Nahuel Freitas and Juan Pablo Paz. Analytic solution for heat flow through a general harmonic network. Physical Review E, 90(4):042128, 2014.

## Appendix A Model and analytic solution

We consider the following Hamiltonian,

 H=HS+HE+Hint, (10)

where is the Hamiltonian of an arbitrary network of quantum harmonic oscillators, and describes a set of uncoupled bosonic reservoirs. is the interaction between system and reservoirs. Therefore,

 HS=12PTM−1P+12XTV(t)X, (11)

where and are vectors whose components are position and momentum operators satisfying the usual commutation relations, and . The matrices and are symmetric and positive definite, and can be time-dependent. The environmental Hamiltonian is

 HE=∑αHE,αHE,α=Nα∑j=1π2α,j2m+mω2α,j2q2α,j

That is, the environment consists of several independent reservoirs formed in turn by a collection of quantum harmonic oscillators of mass . The operator is the position operator of the -th oscillator in the -th environment, and its associate momentum. We consider a bilinear interaction between system and reservoirs through the position coordinates

 Hint=∑α∑j,kCα,jkXjqα,k, (12)

where are time-independent interaction constants.

### a.1 Equations of motion

We now derive, working in the Heisenberg picture, the equations of motion for all the operators involved in the Hamiltonian of Eq. (10). For the system operators the motion equations are

 ˙X =M−1P (13a) ˙P =−V(t)X−∑αCαqα, (13b)

where and are vectors formed with the position and momentum operators of the -th reservoir, respectively. Similarly, the matrix has as elements the interaction constants . If the system has degrees of freedom and the -th reservoir is formed by oscillators, then the matrix has dimensions .

Turning to a description in the phase space, we define the -component vectors and . In terms of , the Eqs. (13a) and (13b) can be written as:

 ˙Z+as(t)Z=∑αCαzα, (14)

where the matrix is defined as

 as(t)=(0−M−1V(t)0), (15)

and the matrix is

 Cα=(00−Cα0) (16)

For the operators corresponding to the -th reservoir the equations of motion in phase space are:

 ˙zα+aα(t)zα=∑α¯CαZ (17)

In this case is the matrix

 aα(t)=(0−1Nα/mmΩ2α0), (18)

where is a diagonal matrix containing the squared frequencies of the oscillators of the -th reservoir. Finally, is:

 ¯Cα=(00−[Cα]T0) (19)

In summary, if we are interested in the dynamics of the system, the following set of linear coupled differential equations must be solved for :

 ⎧⎨⎩˙Z+as(t)Z=∑αCαzα,˙zα+aα(t)zα=¯CαZ (20)

with the initial conditions and .

### a.2 An integro-differential equation for the system

In this section we derive an integro-differential equation describing the dynamics of the system only. We start by considering the Green’s function for the homogeneous equation of motion of the -th reservoir. Such function satisfies:

 ddtgα(t,t′)+aα(t)gα(t,t′)=1Nαδ(t−t′). (21)

With initial conditions , the function encodes the response of the -th reservoir to a delta impulse at time . For this simple case it just represents a rotation in phase space: where is the Heaviside step function and

 e−aαt=(cos(Ωαt)sin(Ωαt)(mΩα)−1−sin(Ωαt)mΩαcos(Ωαt)) (22)

The function is an homogeneous solution of Eq. (17) for . A particular solution is . Therefore, the complete solution of Eq. (17) is

 zα(t)=gα(t,0)zα(0)+∫t0gα(t,t′)¯Cαz(t′)dt′, (23)

which satisfies the required initial condition.

The solution for of Eq. (23) can now be inserted in Eq. (14), the differential equation for the system coordinates. Doing this we obtain:

 ˙Z+as(t)Z−∫t0[∑αCαgα(t,t′)¯Cα]Z(t′)dt′=∑αCαgα(t,0)zα(0) (24)

This is a non-Markovian equation of motion for the system with a source term depending on the operators of the environment at the initial time. Since the equation is linear a general solution can be obtained in terms of the Green’s function of the homogeneous system, as we did before for each reservoir. Before doing that, we define the dissipation kernel.

### a.3 Dissipation kernel

The quantity multiplying in the integrand of Eq. (24) is the dissipation kernel, to which we will refer as , where . It can be written in a more convenient way in terms of the spectral densities of the reservoirs, which are defined below. Explicitly, we have:

 ηα(t,t′)=θ(t−t′)(00ηxxα(t−t′)0) (25)

where the matrix is defined as:

 ηxxα(t)=∫∞0Iα(ω)sin(ωt)dω (26)

The function is the spectral density associated to the -th reservoir. It is defined as follows:

 [Iα(ω)]j,k=Nα∑p=11mω[Cα]jp[Cα]kpδ(ω−ωα,p) (27)

### a.4 Solution of the equation of motion

Using the previously defined dissipation kernel the equation of motion for the system is:

 ˙Z+as(t)Z−∫t0η(t,t′)Z(t′)dt′=∑αCαgα(t,0)zα(0) (28)

We consider the Green function associated with the previous equation. It is such that:

 ∂∂tG(t,t′)+as(t)G(t,t′)−∫t0η(t,τ)G(τ,t′)dτ=12Nδ(t−t′) (29)

with initial conditions . The function is therefore the response of the system to an impulse at time . It fully takes into account the non-Markovian nature of the dynamics and the dissipation induced by the environment. In some cases the function can be computed analytically. In general only a numerical approach is possible. In any case, if is known, the complete solution to Eq. (28) can be obtained. In fact, it is easy to verify that the expression

 Z(t)=G(t,0)Z(0)+∫t0G(t,t′)[∑αCαgα(t′,0)zα(0)]dt′ (30)

is a solution of Eq. (28) and satisfies the required initial condition.

### a.5 Renormalization and damping kernel

It is useful to rewrite the integro-differential equation in Eq. (29) and express the non-Markovian integral term as a functional of the velocity in phase space, , instead of . For that purpose a partial integration must be performed, with the following result:

 ∂∂tG(t,t′)+aR(t)G(t,t′)+∫t0γ(t,τ)∂∂τG(τ,t′)dτ=12Nδ(t−t′) (31)

Note that the matrix , that describes the unitary dynamics of the system, has been renormalized to . The function is such that , and is known as the damping kernel. For it can be calculated as follows:

 γα(t,t′)=(00γxxα(t−t′)0) (32)

where the matrix is defined as:

 γxxα(t)=∫∞0Iα(ω)ωcos(ωt)dω (33)

### a.6 Evolution of the covariance matrix

The covariance matrix of the system at time is defined as

 C(t)=Re[⟨Z(t)Z(t)T⟩]−⟨Z(t)⟩⟨Z(t)T⟩ (34)

where , and is the initial state of the system and reservoirs.We will consider initial states such that and therefore, according to Eq. (30), for all . Inserting Eq. (30) in Eq. (34) the following expression is obtained:

 C(t)=G(t,0)C(0)G(t,0)T+G(t,0)Re[⟨Z(0)β(t)T⟩]+Re[⟨β(t)Z(0)T⟩]G(t,0)T+Re[⟨β(t)β(t)T⟩] (35)

where is the integral term of Eq. (30),

 β(t) =∑αβα(t)zα(0) (36a) βα(t) =∫t0G(t,t′)Cαgα(t′,0)dt′ (36b)

The first term in Eq. (35) is the deterministic propagation of the initial covariance matrix given by the phase space flow . The second and third terms are the propagation of the initial correlations between system and reservoirs. The last term correspond to the noise and diffusion induced by the environment on the system, and in a stable system will dominate the long term behaviour. If there are no system-reservoirs correlations in the initial state, i.e, if for all , then the second and third terms of Eq. (35) vanish for all . In that case the evolution of the covariance matrix is just:

 C(t)=G(t,0)C(0)G(t,0)T+Re[⟨β(t)β(t)T⟩] (37)

Using Eq. (36a) we find the following expression for the diffusive term of the covariance matrix:

 Re[⟨β(t)β(t)T⟩]=∫t0∫t0G(t,t1)⎡⎣∑α,βCαgα(t1,0)Re[⟨zα(0)zβ(0)T⟩]gβ(t2,0)TCTβ⎤⎦G(t,t2)Tdt1dt2 (38)

Now we introduce the condition that in the initial state the reservoirs are in thermal states and uncorrelated with each other. In that case:

 Re[⟨zα(0)zβ(0)T⟩]=δα,βℏ2⎛⎜⎝(mΩα)−1coth(ℏΩα2kBTα)00(mΩα)coth(ℏΩα2kBTα)⎞⎟⎠, (39)

Where is the temperature of the -th reservoir and is the Boltzmann constant. Introducing Eq. (39) in Eq. (38) the following final expression is obtained:

 Re[⟨β(t)β(t)T⟩]=ℏ2∫t0∫t0G(t,t1)ν(t1−t2)G(t,t2)Tdt1dt2 (40)

The matrix function is the noise kernel, with

 να(t)=(000νxxα(t)) (41)

where:

 νxxα(t)=∫∞0Iα(ω)cos(ωt)coth(ℏω2kBTα)dω (42)

### a.7 Asymptotic state for stable systems

In this section we characterize the asymptotic state for driven systems that are exponentially stable, i.e., systems with a Green’s function decaying exponentially with . The previous condition is not always fulfilled for driven systems, even in the presence of strong dissipation, since it is possible, for example, to induce a divergent dynamics by the phenomenon of parametric resonance.

From Eq. (37) it is clear that for stable systems the asymptotic covariance matrix is:

 (43)

since for large . Also, Eq. (40) is equivalent to the following expressions for each block of the asymptotic covariance matrix:

 σxx(t) =ℏ2∫t0∫t0g(t,t1)M−1νxx(t1−t2)M−1g(t,t2)Tdt1dt2 (44a) σxp(t) =ℏ2∫t0∫t0g(t,t1)M−1νxx(t1−t2)M−1∂∂tg(t,t2)TMdt1dt2 (44b) σpp(t) =ℏ2∫t0∫t0M∂∂tg(t,t1)M−1νxx(t1−t2)M−1∂∂tg(t,t2)TMd