Prethermalization of quantum systems interacting with non-equilibrium environments

The usual paradigm of open quantum systems falls short when the environment is actually coupled to additional fields or components that drive it out of equilibrium. Here we explore the simplest such scenario, by considering a two level system coupled to a first thermal reservoir that in turn couples to a second thermal bath at a different temperature. We derive a master equation description for the system and show that, in this situation, the dynamics can be especially rich. In particular, we observe prethermalization, a transitory phenomenon in which the system initially approaches thermal equilibrium with respect to the first reservoir, but after a longer time converges to the thermal state dictated by the temperature of the second environment. Using analytical arguments and numerical simulations, we analyze the occurrence of this phenomenon, and how it depends on temperatures and coupling strengths. The phenomenology gets even richer if the system is placed between two such non-equilibrium environments. In this case, the energy current through the system may exhibit transient features and even switch direction, before the system eventually reaches a non-equilibrium steady state.


I. INTRODUCTION
The standard theory of open quantum systems (OQS) typically considers that the system is coupled to a single reservoir in equilibrium to analyse properties such as decoherence, dissipation and non-Markovianity [1][2][3][4]. A richer situation emerges in the frame of quantum thermodynamics and thermal machines, in which the system is coupled to two or more reservoirs, each of them equilibrated at a different temperature and/or chemical potential [5][6][7][8]. Once the coupling is activated, the open system evolves towards a non-equilibrium steady state that may contain persistent heat, particle or spin currents. An even more involved scenario occurs when the system is coupled to one or more reservoirs that are each of them out of equilibrium and therefore evolve in time. As a consequence, the action of the environment into the system dynamics is no longer encoded in a correlation function that is time-translational invariant, such that α(t, τ ) = α(t − τ ), but rather on a correlation function that depends on both the current time of evolution t and the past times τ .
The motivation to analyze complex environments beyond the standard OQS paradigm of single and multiple equilibrium reservoirs is strong. From an application perspective, out of equilibrium environments that present a temperature gradient can be encountered in electron transfer processes in quantum chemistry and biology [9], in cellular media [10] and even in the thermosynthesis processes that use the solar energy to create chemical compounds [11], to name just a few examples. These types of environments may indeed be driven by an external source, corresponding to other molecular or biological structures or even to the electromagnetic field. Non-equilibrium environments are also present in quantum technological devices, where the quantum system of interest may be directly coupled to an environment that is itself coupled to a second reservoir, thermalized at a different temperature. Such temperature gradient of the different components and subsystems surrounding the quantum system of interest is particularly present in quantum computers [12,13]. Superconducting qubits, for instance, are cooled down to cryostatic temperatures, while their surrounding components, including amplifiers and processing units, as well as the cables and waveguides that connect them to each other and to the qubit, are at higher temperatures the further they are from the circuit.
Describing these situations is of fundamental and timely interest, but it also represents a significant challenge, as the effects of indirect reservoirs on the OQS dynamics can not be captured with a simple Markovian approximation. To this aim, one possibility is to compute the full dynamics, including the system and the environments, and then trace out the environmental degrees of freedom to obtain the OQS reduced dynamics. However, the dimension of the full Hilbert space grows exponentially fast with the number of degrees of freedom, and further, the relevant states may be largely entangled as well, which makes inefficient a direct use of state of the art numerical methods like Monte Carlo [14][15][16] and matrix product states [17][18][19][20][21]. As an alternative, this paper proposes to extend the standard tools of the OQS theory, namely the weak coupling approximation and the master equation approach, to consistently tackle the problem in at least a limit of interest.
To be specific, we consider a two-level quantum system coupled to a first reservoir (RI) that is in turn coupled to a second reservoir (RII). Initially, each reservoir is in a thermal state at a different temperature, respectively T I and T II . We additionally consider that RII induces a Markovian evolution on the modes of RI so that they thermalize efficiently. Therefore, even if RI is initially in thermal equilibrium, the coupling to a second reservoir at a different temperature will drive it away from it, and enforce its evolution towards a new equilibrium state with respect to RII. Thus, the dissipation of the open system will display very rich features reflecting the interplay between two different time scales: thermalization of the system at a temperature T I , and the thermalization to its final equilibrium state with T II . If the conditions of the environment are suitable, and these two time scales are temporally separated, prethermalization [22] of the OQS is observed, which is a stage in which the system remains thermalized at T I .
The plan of the paper is the following: We present the details of our model in Sec. II, while in Sec. III we discuss the master equation that is used to describe the reduced dynamics of the open system. This master equation depends on a set of correlation functions that encode the effects of both reservoirs in the open system, and which are discussed in Sec. IV. Sections V and VI describe the effects of prethermalization when considering a single and two out of equilibrium reservoirs, respectively. Finally, we draw some conclusions in Sec. VII.

II. MODEL WITH TWO INTERACTING ENVIRONMENTS
As is standard in the theory of OQS [3,23,24] we consider that the total evolution of system plus the environment is unitary and described by the Hamiltonian, where H S and H E are the free Hamiltonians of the system and environment, respectively, and H int is the interaction Hamiltonian between system and environment. We model the system as a two level system with the free Hamiltonian where ω 0 is the energy 1 gap between levels. We model the environment to which the system is coupled as an set of open harmonic oscillators that is a first reservoir (RI) of harmonic oscillators where each mode in RI is coupled to a independent reservoir, included in the second reservoir (RII). The Hamiltonian describing this environment is where 1 Throughout this article we consider natural units in which the reduced Plank constant and the Boltzmann constant = k B = 1.  are the free Hamiltonians of RI and RII, respectively, with operators that obey the commutation relations and whose interaction conserves the boson number. The coupling strength between the λ-th oscillator in RI and the k-th oscillator in RII isg λ,k . The system is in a star configuration, i.e. the OQS is coupled to all the λ bosonic operators of RI, and in turn each of these is coupled to a reservoir of harmonic oscillators that is a part of RII, as depicted in Fig. 1. Only the first reservoir couples directly to the OQS, with the interaction Hamiltonian which only considers interactions that conserve the particle number. We take as initial state a tensor product, The initial state of the system can be arbitrary, while the initial states of the reservoirs are assumed to be thermal, possibly at different temperatures, where Z(β i ) = Tr{e −βiH Ri } is the partition function and β i = 1/T i is the inverse temperature of each reservoir, i = {I, II}.
Each reservoir may have a different spectral function depending on their microscopic properties and the problem considered. In our analysis, we consider the Caldeira-Leggett phenomenological model of spectral functions [25], which reads where g i is the strength of the coupling, s i is a factor that takes different values depending on the particular environment that needs to be modelled, and ω ci determines a smooth frequency cut off for the modes of the reservoir.

III. STUDY OF THE SYSTEM EVOLUTION
To obtain a closed equation for the dynamics of the open quantum system, we consider that it is weakly coupled to its environment, which makes it evolve slowly. Thus we can derive a second order weak coupling master equation (ME) for the reduced density matrix of the open quantum system. The derivation of the ME is standard and can be found in numerous works [23,26], where V t O = e iH S t Oe −iH S t is the free evolution of the operator O = {σ + , σ − }, and the correlation functions are defined by with B(t) = e iH E t Be −iH E t the free evolution of the environment operator B = λ g λ a λ . Notice that this equation is second order in the interaction operator B, and that no first order term is present, since it is proportional to Tr E {H int (t)ρ E (0)}, which is null for the initial state defined in Eqs. (8,9). This equation is a time-local ME, since its evolution can be recast in the forṁ where Λ t is a linear map, such that Λ t [ρ(t)] is Hermitian and traceless for any ρ. To fully describe the OQS through the differential equation (11), the correlation functions (12) have to be computed for the initial states ρ E (0) defined in Eqs. (8,9). The following section is devoted to this derivation, but first we rewrite the ME in Eq. (11) under its canonical form.
A. Canonical Form of the ME Any time-local ME equation of the form (13) can be recast into a canonical ME [27], of the form where γ k (t) are the canonical decay rates corresponding to the canonical decoherence channels L k (t), with k = 1, . . . , d 2 − 1, and d the dimension of the Hilbert space of the OQS. H(t) is, in general, not identical to the free Hamiltonian of the system, since the interaction with the environment modifies it. The most common effect is a shift of the natural frequency of the OQS, the so-called Lamb shift. The equation is often written in a more compact form as where the first term represents the unitary evolution of the OQS. The second term in (15) encompasses the dissipative part of the evolution.
Recasting the time-local ME in this form allows us to easily evaluate whether, despite being an approximated equation, it still preserves complete positivity of the evolution. In detail, if the decay rates γ k (t) are non-negative we can ensure that this is the case and that the dynamical map of the OQS is Markovian [27]. The canonical decay rates, and the Lamb shift for our model, are discussed in the next section and in Appendix A.

IV. OUT-OF-EQUILIBRIUM CORRELATION FUNCTIONS AND DECAY RATES
Obtaining the correlation functions (12) requires to compute the time evolution of a λ (t) in the operator B(t) = λ g λ a λ (t). We can simplify this calculation by assuming a large separation of time scales between the second and the first reservoir. Specifically, we consider that the modes of the first reservoir, a λ (t) slowly evolve towards an equilibrium state with respect to the second reservoir, and that this evolution is well described with the Markov approximation. This is discussed in Appendix B, while the computation of the correlation functions is treated in Appendix C. Thus, the correla-tion functions are given by where J i (ω), with i = {I, II}, are the spectral functions of each reservoir, which have the general form (10), and which is proportional to a Lorentzian kernel of width J II (ω)/2, and the function Notice that, even though we can consider that the open system is weakly coupled to RI, and thus its master equation is obtained within a second order perturbation theory, a Markov approximation can not be taken in a straightforward way. The reason is that the correlation functions (16) are no longer dependent on the time difference t − τ , but on both times t and τ such that one can not simply extend the integration limits in Eq. (11) by assuming that the integral kernel decays much faster than the system evolution time-scale, as it is done in the Markov approximation. We observe that the second term of the correlation functions contains the resonant term K(ω, ω ) (see Eq. (17)) with a width proportional to the coupling strength between environments, and centered at ω = ω . Approximating this term by a delta function is consistent with the weak coupling approximation already considered between RI and RII. Using this approach, we obtain an analytical approximation for the canonical decay rates γ ± (t), which correspond to the decoherence channels L ± = σ ± (see appendix A), and which can be split into two contributions, γ ± (t) = γ ST where the terms are labelled in reference to their short time (ST) or long time (LT) dominance. The ST terms are and the LT terms read The validity of approximating Eq. (17) by a delta function is discussed in Appendix D. These decay rates present a very suggestive form: at short times, the LT terms of each decay rate is negligible, while at later times it dominates (see Appendix D for a visual reference). The strength of the decay rates is governed by the spectral function of the first environment, while the second environment spectral function is responsible for the time scales at which each term dominates. With this approximate expression for the decay rates it is possible to prove analytically that indeed the OQS evolves, at long times, to a thermal state at the inverse temperature of the second reservoir β II (see appendix E). Furthermore, since they are non-negative at all times, we can ensure that the ME preserves complete positivity.

V. PRETHERMALIZATION
The decay rates obtained in the previous section already suggest that the evolution of the OQS may exhibit a separation of timescales, such that there may be a transitory state in which the OQS remains close to a thermal state corresponding to the initial temperature of RI, but after some longer time it finally relaxes to a thermal state with the temperature of RII.
This transient effect is an instance of prethermalization, a phenomenon in which the system, after a short time, seems to relax to a state different from the true thermal equilibrium, which is eventually reached after a much longer timescale [22,[28][29][30][31]. The most studied scenario of prethermalization concerns weakly nonintegrable systems, in which an eigenstate of an integrable model is evolved under a quenched Hamiltonian that weakly breaks integrability. The short time dynamics is still determined by almost conserved quantities, and the system arrives to a prethermalized state, but at long times the breaking of the integrability dominates and the system finally thermalizes [32]. The phenomenon has also been studied in the context of OQS in [33,34] and observed experimentally in ultra-cold bosonic atoms [35][36][37].
In our setup, a small coupling to the second reservoir (g II = 0) can play a similar role to the integrability breaking, as it perturbs the thermal equilibrium of the enviroment RI (which would otherwise remain stable). In this way, the initial temperature of RI may determine the short time evolution and the arrival to a prethermal state, while the final equilibrium is determined by RII. We will thus consider that prethermalization has occurred when the system reaches a state, independent of its initial conditions, close to the thermal equilibrium at β I , and this state is mantained for a finite time, before the evolution definitely drives the system to the equilibrium with RII.  29), represented by a red dot. The OQS stays close to this thermal state for some time: The prethermalization time tpr. Afterwards, the OQS starts the evolution towards the thermal state at βII. This corresponds to the displacement of a point (marked as a cross) from the red dot to the green one (representing the thermal state at βII). The arrows connect the points of the upper population with the snapshots of the evolution of the ball of accessible states. The environments parameters are gI = 10 −2 , gII = 10 −5 , sI = sII = 1, ωcI = ωcII = 10, βI = 1 and βII = 0.1 and the system frequency is ω0 = 1.
In order to verify the occurrence of the effect, we analyse the evolution of all possible initial states. We conveniently express the density matrix in terms of the polarization vector, ρ(t) = (I + p(t) · σ)/2 and integrate the time evolution equations (see Appendix A). The formal solution for the polarization vector is where is a rotation matrix, that performs a rotation about the z axis with angular frequencyΩ(t) where Ω(t) is the shifted frequency of the OQS due to the action of the environment (see Appendix A), r(t) = e −Γ(t) is a scaling factor that affects equally all components, with Γ(t) = t 0 dt (γ + (t ) + γ − (t )), where γ + (t) and γ − (t) are the canonical decay rates, and d(t) = (0, 0, c(t)) is a displacement vector in the z direction, with From this result, it is apparent that the effect of the dynamical map on any state is to rotate the polarization vector around the z axis, rescale it by r(t) and add a displacement c(t) along the vertical direction. These transformations are independent of the initial state, hence the space of accessible states, initially described by the volume limited by the Bloch sphere, is isotropically contracted and shifted, and can be characterized by its time dependent radius and center.
We would like to emphasize that the Lamb shift does not play any role in the evolution of the diagonal elements of the reduced density matrix, which ultimately means that it does not affect either the long time thermalization or the prethermalization dynamics. It is encoded in the angular frequency of Eq. (22) and thus it has the effect of rotating the ball of accessible states with an angular velocity different from ω 0 , but does not affect the rescaling and displacement of the whole space.
This representation allows us to understand how fast the memory of the initial state is lost, and in which state the OQS is. For the approximate decay rates Eqs. (19,20) we obtain the following expression for the radius of the ball of accessible states which in the limit J II (ω 0 )t → 0 becomes that is the expression that we would obtain if only RI was considered. This means that the rate at which the volume of the accessible states reduces is mainly governed by RI. A smaller coupling between OQS and RI would cause a slower reduction of the accessible states space. The center of the ball of accessible states is given by the evolved polarization vector of the maximally mixed state, namely the origin of the Bloch sphere, and is thus at c(t) = (0, 0, c(t)) with This expression has no analytic solution, but can be solved in the short time limit (ST), J II (ω 0 )t 1. If the exponentials inside the second factor are Taylor expanded in terms of J II (ω 0 )t and J II (ω 0 )t ≤ J II (ω 0 )t up to first order, the resulting integral is solvable and yields Within this regime, we distinguish two limiting cases which at time t = 0 corresponds to the center of the Bloch sphere.
• When J I (ω 0 )(2n I (ω 0 ) + 1)t 1, the exponential in Eq. (27) vanishes, and this expression becomes such that (0, 0, c ST ) corresponds to the thermal state ρ th S (β I ). This expression holds when in which case the ball of accessible states is centred around the point corresponding to the thermal state of the OQS at β I as long as J II (ω 0 )t 1. Moreover, in this limit the radius of the ball of accessible states Eqs. (24,25) is close to 0, meaning that the state of the OQS is independent of the initial condition and close to the state ρ th S (β I ), which shows that the system thermalizes to β I .
Eq. (30), shows that the condition for the OQS to prethermalize to β I , depends on the relationship of this temperature and the coupling strengths, but is independent of β II . In the next section we analyse how β II affects the prethermalization.  Fig. 2, all initial condition directly tend to the thermal state at βII, which is the asymptotic state of the system. The ball of accessible states contracts around the point representing the thermal state at βII, and after becoming punctual it stays there. There is no prethermalization phenomenon in this case, where the spectral functions parameters are as in Fig. 2 with the exception that gII = 10 −2 .
The long time (LT) limit (J II (ω 0 )t → ∞) of Eq. (26), studied analytically in Appendix F, yields where the point (0, 0, c LT ) corresponds to the thermal state ρ th S (β II ) as the asymptotic state. This asymptotic state was also checked analytically using the approximate decay rates in Appendix E.
To illustrate the above discussion, we display in Figs. 2 and 3 the evolution of the system in two different scenarios. In both cases, the time dependence of the ρ ++ (t) = +| ρ S (t) |+ component 2 of the state of the system is shown for several initial pure states, which allows us to visualize the evolution of the ball of accessible states. In the first case, for β I = 1, β II = 0.1 and g II = 10 −5 we observe prethermalization (Fig. 2), but when the coupling is increased to g II = 10 −2 (Fig. 3), the phenomenon does not appear.
Following our previous considerations, we identify two relevant time scales that govern the OQS evolution in the prethermalization regime of Fig. 2. First, the time t I after which the OQS has evolved to the thermal state at β I . At this time, the space of accessible states has already contracted to a point, so that the state reached is independent of the initial condition. The second time scale t II determines the time required for thermalization to the asymptotic state ρ th S (β II ). If t I is sufficiently smaller than t II , as in Fig. 2, the system first evolves to ρ th S (β I ) (red dot), and stays close to it for a certain time t pr , which we call prethermalization time. After this time, it smoothly evolves to ρ th S (β II ) (green dot). As shown in Fig. 3, when  i.e., T (ρ th S (βI), ρ th S (β i II )), so that when the OQS approaches the asymptotic state, solid lines tend to the dashed lines. Initially the OQS is in a thermal state at βI and departs from it as it evolves. We observe that this departure happens faster for a larger separation between the thermal states of the environments. The red line represents the distance dpr = 10 −2 and the rest of the parameters are as in Fig. 4. the conditions of the problem do not allow for prethermalization, we observe the thermalization of any initial condition directly to the state ρ th S (β II ), without any transitory approach to ρ th S (β I ).

A. Prethermalization Time
To give a more quantitative estimation of the time during which the OQS remains approximately thermalized at the temperature β I , i.e. the prethermalization time t pr , we make use of the trace distance between the evolved state of the OQS, with initial condition ρ S (0) = ρ th S (β I ), and the thermal state ρ th S (β I ). We define t pr as the time elapsed between the time at which the radius of the ball of accessible states has reduced below 10%, and the time at which the above trace becomes bigger than a fixed trace distance d pr . This represents a threshold distance below which two states could not be distinguished. If the order in which these events happen is the opposite, it means that no prethermalization is present.
We can visualize this by looking at the dynamics of the polarization vector corresponding to the density matrix of the system, starting from ρ th (β I ). If that point has been significantly displaced before the ball of accessible states has contracted, then no prethermalization is present: See Figs. 2 and 3 for a visual reference of this criterion. If the trace distance between the thermal state of the system at β I and β II is smaller than d pr , the prethermalization time is not defined, as these two states would not be distinguishable.
With this definition we studied how t pr varies as a function of the initial temperatures of both reservoirs, as well as for different values of the coupling strength between them, i.e. g II . In Fig. 4 we show the prethermalization time as a function of the trace distance for fixed β I varying β II . We observe the prethermalization time to be longer, the closer the two states are. The same can be appreciated in Fig. 5, which shows the calculation of t pr for different separations of the thermal states at β I and β II : When they are closer (orange line) t pr is higher and when they are further apart (blue and green lines) t pr decreases.
We realized that the prethermalization scale for fixed temperatures of both reservoirs depends inversely with g II , i.e., t pr ∼ g −1 II , for small values of g II . This can be clearly seen, for small coupling strenghts, in Fig. 6. In this figure we checked that condition (30) is fulfilled for all points. It was pointed out in [28] that the thermalization time t II scales as g −1 II similarly as t pr , and we checked  that this is true for our case 3 as well. From Figs. 4 and 6 we also observe that the time the OQS stays in the prethermal state is shorter for a higher temperature of RII (lower β II ), while the duration of the prethermalized state increases with β II . A similar qualitative behaviour is observed for varying β I (see right panel of Fig. 6), but the overall prethermalization scale is smaller in case of a larger β I .

VI. COMPOSITE NON-EQUILIBRIUM ENVIRONMENTS
In the scenario discussed in the previous sections, the OQS is expected to reach a steady state at thermal equilibrium, consistent with the temperature of the largest reservoir RII. However, using the same ME formalism it is also possible to construct more complex scenarios in which the steady state of the open system is out of equilibrium, for instance when the system is coupled to two independent environments, as depicted in the left column of Fig. 7. The additional environment gives rise to a new dissipative term in the master equation (15) and a new term in the environment corrected Hamiltonian H(t). The corresponding rates and Lamb shift corrections are computed identically as before. The heat flow between environments produces a change in the energy of the system E S = Tr{H S ρ S (t)}, where H S is the system Hamiltonian. The evolution of this quantity can be expressed by the canonical ME (15) as where the superindex (ν) = {L, R} refers to the left and right environments andH(t) = H S + ν 1 2 ∆ω (ν) (t)σ z . SinceH(t) and H S are both proportional to σ z , the first term vanishes, and the second one defines the heat fluxes from environment (ν) to the OQS. By convention, we consider the heat flux from the right reservoir to be positive, and the one from the left to be negative. Thus, a positive total heat flux indicates a flow of energy from right to left, and vice versa. In the following subsections we first analyze the heat fluxes when the left and right reservoirs are each in a thermal equilibrium state, and then when each of them are out of equilibrium.

A. Heat flux between environments in equilibrium
In the case of equilibrium environments of Fig. 7a, which we depict as single reservoirs on each side of the OQS, any initial state reaches a non equilibrium steady state that depends on the initial state of both environments under our model assumptions. In Fig. 7c we plot the total heat flux J L (t) + J R (t) calculated using Eq. (34) with γ I (ω 0 ) + 1) the decay rates of the spin boson model with one reservoir. We observe that initially, the heat flux depends on the initial condition, but after some time it always converges to the value in the steady state, where ρ qs ++ (β L I , β R I ) is given in Appendix G. When the OQS reaches the asymptotic state there is a constant heat flux from the environment with the higher temperature. In Fig. 7c, which corresponds to β R I < β L I , this is observed by a positive steady state flux.

B. Heat flux between environments that are out of equilibrium
We now consider the case where the OQS is coupled to two out of equilibrium reservoirs, as schematically depicted in Fig. 7b. Fig. 7d shows that the heat flux, independently of the initial condition, is dominated by the temperature gradient between β R I and β L I , while the gradient for β (R,L) II becomes relevant at longer times. The time scale in which each gradient is dominant is determined by g (ν) II . Interestingly, we observe that the interplay between these gradients may even produce a change of sign in the current. This is because we have chosen β R I < β L I , but β R II > β R II , such that the quasi-stationary flux is positive (the right RI is hotter than the left RI), while at long times is negative (since the right RII is colder than the left RII).
Moreover, one can tune these gradients and the couplings g (ν) II to be such that there are two changes of sign in the heat current. This is observed in Fig. 7e, where β R I < β L I < β R II < β L II and g R II > g L II , which leads to an initial and final positive flux (β R i < β L i ). But as g R II > g L II there is some time that the quasi-stationary flux is determined by β L I < β R II , such that the flux during that time is negative.
The stationary and quasi-stationary states for the setup Fig. 7b can be explicitly derived, and are shown in Appendix G.

VII. CONCLUSIONS
We have presented a model to describe an OQS which is coupled to a hierarchy of environments at different temperatures, a situation that can be found in complex environments and interfaces that are present in both natural and quantum technological scenarios. To the best of our knowledge this is the first attempt to account for the presence of external reservoirs that are not directly coupled to the OQS. Although these situations are in principle very complex to analyse, we have shown here that, under certain constraints, one can extract a well-behaved master equation that allows such a description in relevant limits.
In detail, we have considered an open system directly coupled to a reservoir RI, at an inverse temperature β I , that is driven out of equilibrium because of its coupling to a second reservoir RII at β II . With the use of weak coupling and Markovian approximations, we have derived a master equation to describe the evolution of the reduced density matrix of the system, by tracing out the evolution of the environment. Even with these approximations, we were able to observe a rich dynamics of the open system, with the existence of a transitory state, called prethermal state, before the final thermalization, which was found to be determined by the larger reservoir solely. We investigated under which conditions prethermalization is present, and concluded that this state is longer lived when the reservoir RI, directly coupled to the OQS, is hotter and RII colder, as well as when the coupling between reservoirs is the smallest possible. We presented a way to characterize prethermalization that is independent of the initial condition of the OQS, through the evolution of the volume of accessible states.
We have also shown that non-trivial dynamics and competing time scales are also present when we consider two out of equilibrium environments coupled to the system. It is well-known that, in the standard situation where the environments are in equilibrium, a heat flux with a given direction (from the hot to the cold reservoir) is established and prevails at long times. Interestingly, when considering out of equilibrium environments we observe that the time scales induced by different environments may induce that the heat flux switches direction, even more than once.
As shown, the OQS dynamics and its currents do not evolve according to a single time scale, but present a richer dynamics that may be evident in experiments and quantum information processes, particularly at long times. The presence of a prethermalization transitory may be harnessed in quantum technological applications, for instance by considering the initialization protocols of a qubit based on coupling it to a reservoir [38,39]. The added reservoir can potentially be controlled by a second one, according to our scheme, in order to optimize further the protocol. In other words, our work describes the possibility of manipulating and controlling an open system by externally modifying and controlling the reservoir to which it is directly coupled.
Our scheme can be adapted to include more external reservoirs at different temperatures. Multiple layer environments can be found, for instance, in superconducting quantum computers, where qubits are affected not only by surrounding layers cryogenically cooled, but also by outer layers at increasingly higher temperatures. Considering this reservoir structure would allow us to find additional transitory and steady states of the OQS, which can potentially be harnessed and controlled. An interesting subject for further investigation would also be the consideration of the dynamics beyond the weak-coupling approximation, and the inclusion of non-Markovian effects. For the interaction Hamiltonian considered in Eq. (7), the canonical decay rates and decoherence channels of the master equation (11) are where we defined and The operator H(t) is a modification of the free Hamiltonian of the OQS which, in this case, represents a shift of the natural frequency of the system, given by Therefore, this Hamiltonian can be rewritten as where Ω(t) = ω 0 + ∆ω(t), is the shifted frequency of the OQS due to the action of the environment. The ME for the different matrix elements of the reduced density matrix readṡ where ρ ++ (t) = +|ρ S (t)|+ is the upper population and ρ +− (t) = +|ρ S (t)|− is the coherence, in the |± eigenbasis of H S . We made use of the trace preservation of the dynamical map.
where the change of variable τ = t − t has been performed. The above approximation is, in fact, a Markovian approximation for the interaction with RII, since the operatorã λ (t) only depends on t, so that we have neglected its past evolution. This integral can be solved via the Sokhotski-Plemelj theorem by rewriting the second term in the r.h.s. of Eq. (B3) as γ λãλ (t), where we defined the damping constant γ λ = π k |g λk | 2 δ(ω λ,k −ω λ )−i k |g λk | 2 P 1 ω λ,k − ω λ .
(B5) This approximation allows for an exact solution of Eq. (B3), which after undoing the change of variable introduced above leads to a λ (t) = a λ (0)e −(iω λ +γ λ )t + decay rates, depends on the parameters of the spectral functions (10) for both environments. We have numerically checked the accuracy of the approximation for a broad range of parameters that is relevant for our study. This is illustrated in Fig. 8, where we compare both terms of γ + (t) with the exact result (numerical integration of Eq. (A2)). The ST term does not reproduce the oscillations of the exact solution, while the LT term perfectly matches the exact result.
We also considered a more accurate approximation, which consists in taking J II (ω) independent of ω, instead of approximating K(ω, ω 0 ) by a delta function. This allowed to resolve the oscillatory nature of γ ST ± , which is of frequency ω 0 , but the result is not as intuitive as the decay rates (19,20). which in the steady state limit allows us to obtain ρ ss ++ = lim t→∞ γ L + (t) + γ R + (t) γ L + (t) + γ R − (t) + γ R + (t) + γ R − (t) . (G2) By taking the limit, one obtains the following matrix element of the non-equilibrium steady state ρ ss ++ = J L I (ω 0 )n L II (ω 0 ) + J R I (ω 0 )n R II (ω 0 ) J L I (ω 0 )(2n L II (ω 0 ) + 1) + J R I (ω 0 )(2n R II (ω 0 ) + 1) , (G3) while the coherence matrix element in the steady state is null. The rest of the matrix elements can be obtained from the trace preserving property of the evolution. Even though this steady state does not directly represent a thermal state, one can obtain an effective temperature for the OQS, since it is a diagonal state in the basis of H S . The effective temperature of the OQS in the steady state is defined as which is a function of both temperatures β II of the environments. A very similar expression is obtained for the quasi-stationary states ρ qs where the indices i, j refer to whether the corresponding environment is in the prethermal state (i = I) or in the asymptotic state (i = II). This allows us to obtain the quasi-stationary state of the OQS when both environments are prethermalizing (i = j = I), or the state when one has thermalized while the other remains in the prethemalization stage (i = I, j = II, or interchanged). When i = j = II it corresponds to the asymptotic state of Eq. (G3).