Quantum Limit for Driven Linear Non-Markovian Open-Quantum-Systems

The interplay between non-Markovian dynamics and driving fields in the survival of entanglement between two non-degenerate oscillators is considered here. Based on exact analytical results for the non-Markovian dynamics of two parametrically coupled non-degenerate oscillators in contact to non-identical independent thermal baths, the out-of-equilibrium quantum limit derived in [Phys. Rev. Lett. 105, 180501 (2010)] is generalized to the non-Markovian regime. Specifically, it is shown that non-Markovian dynamics, when compared to the Markovian case, allow for the survival of stationary entanglement at higher temperatures, with larger coupling strength to the baths and at smaller driving rates. The effect of the asymmetry of the (i) coupled oscillators, (ii) coupling strength to the baths at equal temperature and (iii) temperature at equal coupling strength is discussed.


Introduction
The survival of quantum features in hot environments is restricted by decoherence [1]. For quantum features such as entanglement to survive in the presence of the environment, the typical energy scale of the system ω must be larger than the energy scale k B T associated to thermal fluctuations [2]. For typical nano-mechanical resonators with oscillation frequencies between 1 MHz and 1 GHz, the quantum nature of these systems is revealed only if the temperature is well below 10 µK and 10 mK, respectively, however, current cryostats operate at best at 15 mK [3].
Recently, it was shown that the presence of driving forces relaxes the traditional criterium ω/k B T > 1 discussed above [4]. To discuss the way how this condition is relaxed, let γ TB be the time scale associated to the non-unitary dynamics induced by the thermal bath and γ p the pumping rate of the driving field. In terms of these time scales, the quantum limit for driven out-of-equilibrium quantum systems reads ω/k B T eff > 1, where T eff = T γ TB /γ p is an effective temperature [4,5]. For pumping rates larger than the decay rates, γ TB /γ p < 1, the effective temperature of the system turns out to be smaller than the thermal equilibrium temperature T . This allows for the survival of entanglement at higher equilibrium temperatures T . For two frequency-degenerate nano-mechanical resonators of frequency ω = 2π × 15MHz and mass 10 −17 kg with a coupling amplitude ∼ 10 −3 mω 2 and a decay rate γ TB = 5 × 10 −5 ω, this condition implies that steady entanglement can be observed at tens of mK [4]. However, entanglement might survive in the range of Kelvin if the frequency or the coupling amplitude can be increased by a factor 10 [4].
Very recently, a non-degenerate version of this physical model was implemented using a strongly inter-coupled two doubly clamped beams [6]. The coupling was mediated by an exaggerated overhang between the clamped beams. The configuration presented in Ref. [6] sustains two spectrally closely-spaced vibration modes at 2π×246 and 2π×262 kHz with quality factors of 1300 and 2200, respectively. In combination with piezoelectric transducers, that are incorporated directly into the mechanical elements, it provides the key to realizing efficient parametric down-conversion [6].
Moreover, there is by now experimental evidence that the non-unitary dynamics of micro-resonators is driven by non-Markovian dynamics [7]. This in turn suggests that a similar behaviour may be encountered at the nano-mechanical level, a realm where the low temperature condition already may introduce non-Markovian correlations [8]. Motivated by the experimental results discussed above and by the recent interest and important role of non-Markovian dynamics in, e.g., biological systems [9][10][11][12], quantum metrology [13,14], foundations of quantum thermodynamics [8,15] and nuclear reactions [16], and by the intricate and delicate interplay between non-Markovian dynamics and driven fields in optimal control scenarios [17,18], the quantum limit derived in Ref. [4] is extended here to the case of non-Markovian dynamics.
Specially, the non-Markovian entanglement dynamics of a coupled of non-degenerate oscillators parametrically coupled and in contact to independent non-identical thermal baths are analytically solved by means of the influence functional approach of Feynman and Vernon [19]. Results derived here are valid for any parameter regime and allow for predicting that, when compared with the Markovian case, non-Markovian dynamics (i) decrease the time needed to generate entanglement, (ii) increase the temperature and the coupling-strength-to-the-environment limits at which steady state entanglement can be found, (iii) decrease the pumping rate for reaching a particular amount of entanglement and (iv) relax the resonant driving condition.

A Paradigmatic Model for Several Physical Systems
The model described below is capable of describing a large variety of physical systems such as coupled trapped ions [20], coupled membranes or mechanical oscillators [21]. Specifically, two parametrically coupled harmonic oscillators with different masses m α and frequencies ω α are considered where c(t) is an arbitrary coupling function between the oscillators. This time dependence is critical for the generation and maintenance of entanglement [4,17].
To describe the interaction with their surroundings (dissipative and decoherencing effects) in a rigorously way, the system-bath model [19,[22][23][24][25] in the context of the Caldeira-Leggett model is considered here. To avoid extra correlations between the oscillators from sharing a common thermal bath [26,27], the two oscillators are coupled to independent thermal baths with different power noise [different spectral density J(ω) and different temperature T ], see below. The Hamiltonian of the bathsĤ B , including the interaction with the system of interestĤ I , is then given by where the coefficients c α,k are the coupling constants among the oscillators of the system of interest and each of the modes of their own thermal baths. It can be seen that the interaction is bilinear in the position operators of the systems and the thermal baths. This implies considering only the linear response of the thermal baths to the influence of the system. This consideration is valid only for the case of geometrically macroscopic thermal baths, which leads to a weak interaction among the oscillators in the system and each one of the oscillators comprising the baths [22,24]. Note that in the interaction Hamiltonian (2), there are two terms that only depends on q α and on the coupling constants to the thermal baths. Those terms are included to compensate the renormalization of the harmonic potentials in the system of interest by the presence of the thermal baths [22,24,28]. By considering these terms, it is ensured that the minimum of the global Hamiltonian with respect to the system-of-interest coordinates is determined only by the potentials in the system [28].

Analytic Exact System Dynamics
The dynamics of the coupled oscillators are solved by means of the influence functional approach [19]. Details can be found in Appendix A below. In solving for the dynamics, the initial density operator of the oscillator coupleρ S and their thermal bathsρ TB 1 and ρ TB 2 are assumed to factorize, i.e., that is, it is supposed that at time t = 0 there are no initial correlations between the subsystems of the overall system. However, this assumption may not always valid because in many applications both, the degrees of freedom of the system of interest and the environment to which it is attached form part of the same system. Thus, it is possible that the initial correlations are not available to the experimentalist [25]. Although, initial conditions may be relevant in the generation of new control strategies in open quantum systems [29][30][31][32], they are not considered below for the sake of concreteness. According to the influence functional approach, at time t the matrix elements q 1+ , q 2+ ρ S (t) q 1− , q 2− of the density operatorρ S (t) are given by where q ± = (q 1± , q 2± ) and the propagating function of the reduced density matrix J(q + , q − , t; q + , q − , 0) reads F[q + , q − ] denotes the influence functional [19,33] and is given by [see Appendix A] The dissipation kernel γ α (s) and noise kernel K α (s) are defined in terms of the spectral density J α (ω) as where β α = 1/k B T α , being k B the Boltzmann constant and T α the temperature of each of the thermal baths. The spectral density J α (ω α ) = π N k=1 c 2 α,k 2m α,k ω α,k δ(ω α −ω α,k ) comprises all the information of the bath that is needed to account for its influence on the system.
In deriving an exact closed expression for the propagating function, it is necessary to evaluate the four-fold path integral in Eq. (5). The exact path integration is performed by taking advantage of the linearity of the system and by using standard techniques [22,28]. Details on the derivation can be found in Appendix A. The propagating function is conveniently written as with Q α = 1 2 (q α+ + q α− ), q α = q α+ − q α− and N (t) is a normalization factor that can be determined from the conservation of the normalization of the density matrix, trρ S = 1, the matrix elements A ij are defined in Eq. (A.49). The new coordinates Q α and q α satisfy the following equation of motion with the boundary conditions for Q α (0) = Q α , Q α (t) = Q α , q α (0) = q α , q α (t) = q α . In the original frequency-and-mass degenerate formulation of the driving-assistedhigh-temperature-entanglement scenario [4], the spectral densities J 1 (ω) = J 2 (ω) = J(ω) were taken in the Ohmic form J(ω) = mγω with γ denoting the standard coupling to the environment constant. This leads to local-in-time dissipation in Eqs. (11) provided by the fact that γ 1,2 (s) = 2γδ(s), and induces larger decay rates for the loss of entanglement (see below). To overcome this and to have a more general and complete characterization of the dynamics that allows for a closer description of experimental realization [6,7], the influence of the bath on the system is characterized here by the spectral density where Ω α denotes a finite cut-off frequency. By replacing the last expression in Eq. (7), this spectral density generates γ α (s) = γ α Ω α exp (−Ω α |s|) , which leads to memory effects in the dissipation for times s < τ α = Ω −1 α in the equations of motion (11).
For this spectral density K α (s) = m α γ α Ω 2 α ( β α ) −1 ∞ n=−∞ [Ω α exp (−Ω α |s|) − |ν α,n | exp (−|ν α,n ||s|)] Ω 2 α − ν 2 α,n −1 , being ν α,n = 2πn( β α ) −1 Matsubara's frequencies [22,23,28]. For the degenerate Ohmic situation in Ref. [4], the solution to Eqs. (11) can be expressed in terms of the solutions of Mathieu's oscillator (see also Ref. [34]). The non-Ohmic character of the spectral density in Eq. (12) prevents here the formulation of the solution of Eqs. (11) in a similar fashion. However, their linear character allows for expressing the formal solution in the form where this set of sixteen auxiliary functions {U i , V i , u i , v i } is obtained by numerical integration of the associated set of second order non-local-in-time differential equations that arises for {U i , V i , u i , v i } after plugging (13) into Eq. (11). Because of the timereversed character of the limits in the integral contribution to the equation of motion of coordinate q α in Eq. (11), special care must be exercised in the numerical integration. Note that a direct numerical integration of (11) would not allow for deriving the analytic result for the propagating function in Eq. (9).

Entanglement Quantification and Covariance Matrix Elements
Due to the linearity of the system's Hamiltonian (1) and the Gaussian character of the propagating function of the reduced density matrix in Eq. (9), every initial Gaussian state evolves into another Gaussian state. For the present kind of bipartite system of continuous variables in a Gaussian state, entanglement can be easily quantified in terms of the logarithmic negativity [35]. This quantity gives a characterization of the amount of entanglement which can be distilled into singlets.
To quantify entanglement in this case, Gaussian continuous variable states, only the covariance matrix σ is needed. It reads σ ij = 1 2 ξ iξj +ξ jξi − ξ i ξ j , wherê ξ = (q 1 ,q 2 ,p 1 ,p 2 ) andq 1 ,q 2 ,p 1 ,p 2 are the position and momentum operators of the oscillators in the system of interest characterized by the Hamiltonian (1).

Entanglement Quantification
The logarithmic negativity is defined as where l i 's are the symplectic eigenvalues of the covariance matrix. They are the normal eigenvalues of the matrix −iΣσ, with Σ = 0 I 2 −I 2 0 the symplectic matrix and I 2 the identity matrix of dimension 2. For separable states,ρ S = i p iρ , the logarithmic negativity of the system is zero and each oscillator can be described independently. For continuous variable systems, entanglement has as upper limit the maximally entangled EPR wave-function with E N → ∞. Hence, the amount of entanglement measured by the logarithmic negativity is unbounded from above.

Covariance matrix elements
To calculate any mean value of any operator associated with the observables of one of the oscillators in the system of interest, it is necessary to find the reduced density matrix associated to that oscillator. This is obtained by tracing out the reduced density matrix over the coordinates of the other oscillator, i.e.,ρ S 1,2 = tr S 2,1ρ S . For instance, the matrix elements of the reduced density operator associated with the first oscillator are given by while for the second oscillator To find the first and second moments of each oscillator, it is necessary to perform the following integrals: To find the covariances between the position and/or momentum operators between the oscillators, it is necessary to integrate out For the case when no initial correlations exist between the oscillators in the system of interest, namely,ρ S (0) =ρ S 1 (0) ⊗ρ S 2 (0), the exact analytic expression for the first and second moments as functions of the set of auxiliary functions {U i , V i , u i , v i } can be found at http://gfam.udea.edu.co/ lpachon/scripts/oqsystems. Note that the expressions for the first and second moments (15) and the mixed moments (16) were calculated for arbitrary driving forces c(t). This allows for using expressions (15) and (16) in the context, e.g., of optimal control of sideband cooling of nanomechanical resonators [36].

Initial Gaussian states for the simulations
To obtain specific results, both oscillators are assumed in a general Gaussian state, and therefore, the initial density matrix for each oscillator can be expressed in terms of the coordinates Q α and q α as qp are the first moments of position and momentum, and the variance of the position, momentum and position-momentum, respectively, for the α-th oscillator.

Entanglement Dynamics for Symmetric Thermal Baths
To characterize the influence of the non-Markovian dynamics in reaching a different quantum limit, symmetric thermal bath are considered at this point, i.e., γ 1 = γ 2 = γ, Ω 1 = Ω 2 = Ω and T 1 = T 2 = T . Moreover, to compare with results in Ref. [4], c(t) is chosen here as where ω d denotes the frequency of the driving field and c 1 its constant amplitude.
Although results presented below are particular to the spectral density in Eq. (12), they encompass the most basic feature of non-Markovian dynamics, namely, a non-flat power noise [8].
Under the conditions stablished above, and for degenerate oscillators ω 1 = ω 2 = ω and m 1 = m 2 = m, entanglement is quantify below by means of the logarithmic negativity introduced in Sec. 3.1. The undriven non-Markovian dynamics for this degenerate case was previously analyzed in Ref. [37] and it was found that when ω falls inside the spectral density (non-Markovian case), entanglement persists for a longer time than in the Markovian case. For the present driven case, the dynamic features provided by Markovian [Ω α = 20ω in Eq.   An additional dynamic feature present in Figs. 2 and 1 is that entanglement is generated at shorter times under non-Markovian dynamics than in the Markovian case. Since the rate of the incoherent processes decreases by non-Markovian dynamics (see Sec. 4.4), then it is natural to expect that the driving force needs to preform less work to squeeze the normal modes of the oscillators when the dynamics are non-Markovian. Figure 2 depictes the time dynamics of entanglement for a variety of coupling rates γ at fixed temperature. In analogy to the case discussed in Fig. 1, non-Markovian dynamics are able to support the same amount of steady state entanglement at twice the coupling rate γ than the corresponding Markovian case. Because of the non-Markovian character of the dynamics, simulations over several periods of the driving force are very expensive on computational terms, the amplitude strength c 1 takes a rather large value (c 1 = 0.2mω 2 ) so that the generation of entanglement occurs after a few periods of driving. However, the effects discussed above are clearly present for smaller values of the amplitude strength c 1 (see Sec. 4.4).

Entanglement dynamics as a function of the initial state
One of the most attractive features of the generation of entanglement by driving forces in the presence of non-unitary dynamics is that the system reaches the same amount of stationary entanglement independently of the its initial state [4]. Despite of the non-Markovian dynamics and its associated dependence of the history of the system evolution, this feature remains present here. Figure 3 shows the time evolution of the logarithmic negativity for a variety initial states given by Eq. (17). There, it is clear that all of them reach the same amount of stationary entanglement.

Non-Markovian Quantum Limit
If the characteristic frequency of a given system is denoted by ω, the survival of quantum features such as entanglement can be predicted for the parameter relations satisfying the condition S(ω) < ωµ, being S(ω) = 1 2m J(ω) coth( 1 2 ωβ) the power noise and µ the pumping rate. Specifically, at high temperate ωβ 1, entanglement survives in the steady state if Because of the various time/energy scales involved in this non-trivial out-of-equilibrium situation, the impact of non-Markovian dynamics can understood in manifold ways, specifically, it can be effectively ascribed to each time/energy scale independently. In doing so, the case of the spectral density in Eq. (12) is discussed below and three non-Markovianly scaled parameters, T nM , γ nM and µ nM , are introduced. After plugging the spectral density (12) in the non-Markovian quantum limit given by Eq. (19), a new effective temperature can be defined such that the quantum limit can be cast in the form found in Ref. [4], namely, k B T nM / ω ≤ µ/γ, but with T nM instead of T . It is also possible to define an effective coupling constant For Ω ∼ ω, T nM ∼ 1 2 T and γ nM ∼ 1 2 γ. This scaling of the temperature or the coupling constant explains the results depicted in Figs. 1 and 2, respectively. Alternatively, the non-Markovian scaling factor 1 − 1 1+Ω 2 /ω 2 can be assigned to a third energy scale. Define so that µ nM ≥ µ provided by the fact that 1 − 1 1+Ω 2 /ω 2 ≤ 1. To be more concrete, note that the particular situation analyzed in Ref. [4] corresponds to the case J(ω)/m = γω and µ corresponds to the imaginary part of the associated Mathieu coefficient. Assuming that to leading order in the coupling c 1 [see Eq. (18)] the imaginary part of the Mathieu coefficient can still be expressed by µ ∼ c 1 /4ω, non-Markovian dynamics can be seen as effectively enhancing the coupling between oscillators. This implies that under non-Markovian dynamics, the amplitude of the driving force needed for the entanglement to survive is clearly smaller than in the Markovian case.

Entanglement Dynamics for Asymmetric Thermal Baths
The presence of driving forces above already placed the system into a nontrivial out-ofequilibrium situation. Another nontrivial out-of-equilibrium situation that this system may encounter is the case of environments at different temperatures T 1 = T 2 or resonators with different couplings constants γ 1 = γ 2 . Motivated by the possible role that heat currents may play in the generation entanglement, these two situations are considered below. Fig. 1 it is clear that the lower the temperature, the larger the amount of steady state entanglement that is reached. However, because heat transfer between thermal baths at different temperature is assisted here by the interaction between oscillators, asking for (i) the possible role of heat transfer in preserving/destroying quantum correlations between oscillators and for (ii) the classical/quantum nature of possible correlations established by heat fluxes at high/low temperature among the oscillators, are legitime questions. For concreteness of the present work, these concerns are analyzed in a separated contribution and here interest is restricted to the amount of entanglement that can be reached for different temperature ratios. Specifically, Fig. 4   features in this figure are the strong dependence of the time at which entanglement is generated on the temperature ratio and the dependence of entanglement on the absolute value of the temperature of the baths and not only on the temperature difference. This indeed motivates a comprehensive analysis of entanglement dynamics and heat fluxes.

Thermal baths with different decay rate: γ 1 = γ 2
Because the effective temperature at which each oscillator thermalizes is a function of the power noise of the bath and therefore of the coupling constant [8], another interesting situation from a thermodynamic point of view is the case when the coupling constants are different. Figure 5   entanglement is generated on the coupling ratios, the behaviour is essentially the same as in Fig. 4. Below, the case of non-degenerate oscillators is analyzed and the effect of non-Markovian dynamics in the resonance condition ω d = ω 1 + ω 2 is discussed.

Entanglement Dynamics for Asymmetric Oscillators
Since the case of degenerate oscillators lacks of experimental relevance [38] and the generation of squeezing in mechanical setups was already achieved for non-degenerate oscillators [6], the dynamics of entanglement is considered next for different masses and frequencies. This last situation complements the analysis of the influence of asymmetries started in the previous section with the cases T 1 = T 2 and γ 1 = γ 2 .
The undriven non-Markovian dynamics for the degenerate and non-degenerate cases were previously analyzed in Ref. [38]. In particular, Ref. [38] addresses the effect of the resonance condition for degenerate oscillators and its relationship with the possibility of preserving entanglement at asymptotic times.
analyzed for a variety of frequency ratios ω 2 /ω 1 for fixed ω 1 . To isolate the effect of non-Markovian dynamics, parameters are chosen so that no entanglement is found in the Markovian degenerate case with ω d = 2ω 1 . Figure 6 not only shows that non-Markovian dynamics support the creation and survival of steady entanglement for the degenerate case with "resonant driving", but also over a broad range of frequency detuning (∼ 7%) and with "non-resonant" driving. This feature adds to the known robustness of the squeezing-entanglement generation against small detuning form the resonance frequency.
In the experimental work on the generation of two-mode squeezing reported in Ref. [6], ω 2 ∼ 0.939ω 1 with ω d = ω 1 + ω 2 and γ ∼ 10 −3 , so that the for non-Markovian scenario and pumping rate in Fig. 6, entanglement can be reached at ∼ 37 mK. For quality factors two orders of magnitude larger, entangled modes could be found at 3.7 K. Because entanglement is not found for this set of parameters in the Markovian case, note that reaching these high temperatures is supported by the non-Markovian character of the dynamics. Figure 7 depicts the logarithmic negativity for a variety of mass ratios m 2 /m 1 for fixed m 1 . As it is expected, the smaller the ratio m 2 /m 1 is, the more effective the driving field is in generating entanglement out of the modulation of the coupling strength. Note that the smaller the mass ratio is, not only the larger the value of E N is, but also the shorter the time at which entanglement is generated. This is consistent with the intuitive ideas about effectiveness of the driving field in the presence of lighter masses.

Discussion
For highly symmetric cases, there is evidence that non-Markovian dynamics may allow for the survival of entanglement at temperatures higher than the corresponding Markovian case provided by the interaction with a common bath [27,38]. For purely dephasing baths, non-Markovian dynamics may allow for larger values of entanglement [39]. Based on analytic exact results for the non-Markovian and out-of-equilibrium dynamics of two non-degenerate parametrically coupled harmonic oscillators, it is shown here that in the presence of time-dependent external fields, non-Markovian dynamics support the generation of out-of-equilibrium steady state entanglement at higher temperatures, larger coupling-to-the-environment constants and lower pumping rates than in the Markovian case. This delicate interplay between driving fields and non-Markovian dynamics sets a new quantum limit which incorporates the main time and energy scales of the physical systems at hand.

Appendix A. Influence Functional and Propagating Function of the System
The starting point in the influence-functional theory is considering the density operator of the global system at time t, that in terms of the initial density operatorρ(0) is given byρ whereĤ =Ĥ S +Ĥ B +Ĥ I is the Hamiltonian of the global system withĤ S ,Ĥ B andĤ SB given by Eqs. (1) and (2). Physically, expression (A.1) implies that the global system is considered as a closed one, and therefore, it is possible to evolve it in time by means of an unitary operator. In the coordinate representation, the density operatorρ(t) reads where q α± stands for the coordinates of the oscillators in the system of interest, and q α± = (q α,1± , q α,2± , . . . , q α,N ± ) the coordinates of the thermal baths. The matrix elements for the temporal evolution operators are known as the propagating kernel [19,22,25,28] and is given by This object evolves the overall system "forward in time", that is why the "+" subscript is used. An analogous expression stands for the matrix elements q 1− , q 2− , q 1− , q 2− e i Ĥ t q 1− , q 2− , q 1− , q 2− that evolves the overall system "backward in time" and this is reflected by the "−" subscript in the coordinates. The path integrals in the propagating kernels must be evaluated over the paths q α± and q α± that satisfy the following boundary conditions In equation (A.3), S stands for the action for the global system, defined as usual: ds L S (q 1 (s),q 2 (s), q 1 (s), q 2 (s)) + t 0 ds L B (q 1 (s),q 2 (s), q 1 (s), q 2 (s)) where L S , L B and L I stand for the Lagrangian of the system of interest, the thermal baths and the interaction between these subsystems, respectively. The Lagrangian have the following form By replacing the expression (A.3) into (A.2), the matrix elements of the density operator of the global system read This expression describes the dynamics of the global system. However, all this information is not necessary. The only information that is relevant for the present case is that of the system of interest under the influence on the environment. Therefore, the relevant object here is the reduced density matrix, which is obtained by tracing out over the degrees of freedom of the thermal baths in Eq. (A.9), i.e., At this point, it is assumed that the initial density operator of the global system factorizes,ρ(0) =ρ S (0) ⊗ρ B 1 (0) ⊗ρ B 2 (0). In the position representation,ρ(0) reads By replacing this expression in Eq. (A.10), the reduced density matrix is found to read denotes the propagating function of the reduced density matrix and the object F[q 1+ , q 2+ , q 1− , q 2− ] is the influence functional [19,33] given by (A.14) Appendix A.1. Derivation of the influence functional To derive the influence functional, it is useful to use one of the properties of the influence functional, specfically, If a system is interacting simultaneously with two uncoupled and independents environments A and B, with no initial correlations between them, then [19] Using the above property, it is possible to express the influence functional in (A.14) as the product of two influence functionals, one for each oscillator in the system of interest, namely, Having in mind that the oscillators encompassing each thermal bath are independent among them, the property of the influence functional stated above can be used once more so that each influence functional in Eq. (A.16) is given by where each F k [q α+ , q α− ] describes the influence of each oscillator in one thermal bath on the oscillator in the system of interest. Each of these influence functionals reads Bα (q α,k+ , q α,k− , 0) denotes the initial density matrix for the k th oscillator in the thermal bath. For the case of a thermal bath at thermal equilibrium at a temperature T α , it reads (see, e.g., Ref. [28]) The next step is to evaluate the path integrals in Eq. (A. 19). In doing so, note that the global system under study is linear, and therefore, it is only necessary to evaluate the actions S that can be obtained from the Lagrangians in (A.7) and (A.8). The trick for solving this differential equation consists in treating the system coordinate q α as if it were a given function of time. So that, a differential equation for a driven harmonic oscillator is obtained. Therefore, the path integrals in (A.19) correspond to the kernel for a driven harmonic oscillator (see, e.g., Ref. [28]), except for the term This contribution can be taken out of the path integral because it does not contain any term that depends on the classical paths of the bath oscillators. Having in mind the boundary conditions in Eq. (A.4), taking q α− (t) = q α+ for the tracing operation and using the propagation kernel for a driven harmonic oscillator, the path integrals in (A. 19) are readily given by m α,k ω α,k 2π sinh( β α ω α,k ) m α,k ω α,k 2πi sin(ω α,k t) m α,k ω α,k i 2π sin(ω α,k t) After performing the integrations over q α,k+ , q α,k+ , q α,k− , the influence functional reads The above expression describes the influence of the k − th oscillator in the α − th thermal bath on the α − th oscillator in the system of interest. By replacing this expression in (A.18), the complete expression for the influence functional is obtained for one thermal bath. Specifically, As it is customary in the literature, the limit to the continuum in the spectrum of the bath is performed by means of the spectral density J α (ω α ), which is defined as Therefore, for a discrete set modes in the thermal baths, the above spectral density is made up of Dirac deltas. However, for the set of oscillators in the thermal baths to behave as a formal thermal bath, it is assumed that the frequencies ω α,k are so dense as to form a continuous spectrum. Thus, in the continuous limit, the spectral density can be represented as a continuous and smooth function on the frequency ω α [22]. This allows for defining the dissipation γ α (s) and noise K α (s) kernels in terms of the spectral density as By replacing Eq. (13) into the last expression, the propagating function reads J(Q , q , t; Q , q , 0) = 1 N (t) exp i q 1 U 1 (t, t)Q 1 +U 2 (t, t)Q 1 +U 3 (t, t)Q 2 +U 4 (t, t)Q 2 − q 1 U 1 (t, 0)Q 1 +U 2 (t, t)Q 1 +U 3 (t, 0)Q 2 +U 4 (t, 0)Q 2 + q 2 V 1 (t, t)Q 2 +V 2 (t, t)Q 2 +V 3 (t, t)Q 1 +V 4 (t, t)Q 1 −q 2 V 1 (t, 0)Q 2 +V 2 (t, 0)Q 2 +V 3 (t, 0)Q 1 +V 4 (t, 0)Q 1 × exp − 1 t 0 ds s 0 du K 1 (u − s) [u 1 (t, s)q 1 + u 2 (t, s)q 1 + u 3 (t, s)q 2 +u 4 (t, s)q 2 ] × [u 1 (t, u)q 1 + u 2 (t, u)q 1 + u 3 (t, u)q 2 + u 4 (t, u)q 2 ] respectively. Here φ I and φ R are the imaginary and real components, respectively, of the phase in the propagating function. In terms of the matrix elements A ij and B ij , the propagating function can be written as (A.50)