Quantum limit for driven linear non-Markovian open-quantum-systems

The dynamics of the coupled oscillators are solved by means of the influence functional approach[23]. Details can be found in appendix below. In solving the dynamics, the initial density operator of the coupled oscillators ${{\hat{\rho }}_{{\rm S}}}$ and their thermal baths ${{\hat{\rho }}_{{\rm T}{{{\rm B}}_{1}}}}$ and ${{\hat{\rho }}_{{\rm T}{{{\rm B}}_{2}}}}$ are assumed to factorise; i.e.,

$$\begin{eqnarray}&&\hat{\rho }(0)={{\hat{\rho }}_{{\rm S}}}(0)\;\otimes \;{{\hat{\rho }}_{{\rm T}{{{\rm B}}_{1}}}}(0)\;\otimes \;{{\hat{\rho }}_{{\rm T}{{{\rm B}}_{2}}}}(0);\end{eqnarray} \tag{ 3 }$$

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 be 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 [29]. Although initial conditions may be relevant in the generation of new control strategies in open quantum systems [33–36], they are not considered below for the sake of concreteness.

According to the influence functional approach, at time t the matrix elements $\left\langle q_{1+}^{\prime\prime },q_{2+}^{\prime\prime } \right|{{\hat{\rho }}_{{\rm S}}}(t)\left| q_{1-}^{\prime\prime },q_{2-}^{\prime\prime } \right\rangle$ of the density operator ${{\hat{\rho }}_{{\rm S}}}(t)$ are given by

$$\begin{eqnarray}&&\left\langle {\bf q}_{+}^{\prime\prime } \right|{{\hat{\rho }}_{{\rm S}}}(t)\left| {\bf q}_{-}^{\prime\prime } \right\rangle =\int _{-\infty }^{\infty }{{{\rm d}}^{2}}q_{+}^{\prime }{{{\rm d}}^{2}}q_{-}^{\prime }J({\bf q}_{+}^{\prime\prime },{\bf q}_{-}^{\prime\prime },t;{\bf q}_{+}^{\prime },{\bf q}_{-}^{\prime },0)\left\langle {\bf q}_{+}^{\prime } \right|{{\hat{\rho }}_{{\rm S}}}(t)\left| {\bf q}_{-}^{\prime } \right\rangle ,\end{eqnarray} \tag{ 4 }$$

where ${{{\bf q}}_{\pm }}=({{q}_{1\pm }},{{q}_{2\pm }})$ and the propagating function of the reduced density matrix $J({\bf q}_{+}^{\prime\prime },{\bf q}_{-}^{\prime\prime },t;{\bf q}_{+}^{\prime },{\bf q}_{-}^{\prime },0)$ reads

$$\begin{eqnarray}&&J({\bf q}_{+}^{\prime\prime },{\bf q}_{-}^{\prime\prime },t;{\bf q}_{+}^{\prime },{\bf q}_{-}^{\prime },0)=\int {{\mathcal{D}}^{2}}{{q}_{+}}{{\mathcal{D}}^{2}}{{q}_{-}}{\rm exp} \left\{ \frac{{\rm i}}{\hbar }\left( {{S}_{{\rm S}}}[{{{\bf q}}_{+}}]-{{S}_{{\rm S}}}[{{{\bf q}}_{-}}] \right) \right\}\mathcal{F}[{{{\bf q}}_{+}},{{{\bf q}}_{-}}].\end{eqnarray} \tag{ 5 }$$

$\mathcal{F}[{{{\bf q}}_{+}},{{{\bf q}}_{-}}]$ denotes the influence functional [23, 37] and is given by [see appendix]

$$\begin{eqnarray}\begin{array}{rcl} \mathcal{F}[{{{\bf q}}_{+}},{{{\bf q}}_{-}}] & = & \mathop{\prod }\limits_{\alpha =1}^{2}{\rm exp} \left( -{{\frac{{\rm i}}{\hbar }}^{{}}}\frac{{{m}_{\alpha }}}{2}\left\{ (q_{\alpha +}^{\prime }+q_{\alpha -}^{\prime })\int _{0}^{t}{\rm d}s\;{{\gamma }_{\alpha }}(s)\left[ {{q}_{\alpha +}}(s)-{{q}_{\alpha -}}(s) \right] \right. \right. \\ {} & {} & +\left. \left. \int _{0}^{t}{\rm d}s\;\int _{0}^{s}{\rm d}u\;{{\gamma }_{\alpha }}(s-u)\left[ {{{\dot{q}}}_{\alpha +}}(u)+{{{\dot{q}}}_{\alpha -}}(u) \right]\left[ {{q}_{\alpha +}}(s)-{{q}_{\alpha -}}(u) \right] \right\} \right) \\ {} & {} & \times {\rm exp} \left\{ -\frac{1}{\hbar }\int _{0}^{t}{\rm d}s\;\int _{0}^{s}{\rm d}u\;{{K}_{\alpha }}(u-s)\left[ {{q}_{\alpha +}}(s)-{{q}_{\alpha -}}(s) \right]\left[ {{q}_{\alpha +}}(u)-{{q}_{\alpha -}}(u) \right] \right\}. \\ \end{array}\end{eqnarray} \tag{ 6 }$$

The dissipation kernel ${{\gamma }_{\alpha }}(s)$ and noise kernel ${{K}_{\alpha }}(s)$ are defined in terms of the spectral density ${{J}_{\alpha }}(\omega )$ as

$$\begin{eqnarray}&&{{\gamma }_{\alpha }}(s)=\frac{2}{{{m}_{\alpha }}}\int _{0}^{\infty }\frac{{\rm d}{{\omega }_{\alpha }}\;}{\pi }{{J}_{\alpha }}({{\omega }_{\alpha }}){\rm cos} ({{\omega }_{\alpha }}s),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{{K}_{\alpha }}(s)=\int _{0}^{\infty }\frac{{\rm d}{{\omega }_{\alpha }}\;}{\pi }{{J}_{\alpha }}({{\omega }_{\alpha }}){\rm coth} \left( \frac{\hbar {{\beta }_{\alpha }}{{\omega }_{\alpha }}}{2} \right){\rm cos} ({{\omega }_{\alpha }}s),\end{eqnarray} \tag{ 8 }$$

where ${{\beta }_{\alpha }}=1/{{k}_{{\rm B}}}{{T}_{\alpha }}$, with k_B being the Boltzmann constant and ${{T}_{\alpha }}$ the temperature of each of the thermal baths. The spectral density ${{J}_{\alpha }}({{\omega }_{\alpha }})=\pi \sum _{k=1}^{N}\frac{c_{\alpha ,k}^{2}}{2{{m}_{\alpha ,k}}{{\omega }_{\alpha ,k}}}\delta ({{\omega }_{\alpha }}-{{\omega }_{\alpha ,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 equation (5). The exact path integration is performed by taking advantage of the linearity of the system and by using standard techniques [26, 32]. Details on the derivation can be found in appendix. The propagating function is conveniently written as

$$\begin{eqnarray}\begin{array}{rcl} J({\bf Q}^{\prime\prime} ,{\bf q}^{\prime\prime} ,t;{\bf Q}^{\prime} ,{\bf q}^{\prime} ,0) & = & \frac{1}{N(t)}{\rm exp} \left\{ \frac{{\rm i}}{\hbar }\mathop{\sum }\limits_{\alpha =1}^{2}\left[ q_{\alpha }^{\prime\prime }{{{\dot{Q}}}_{\alpha }}(t)-q_{\alpha }^{\prime }{{{\dot{Q}}}_{\alpha }}(0) \right] \right. \\ {} & {} & \left. -\frac{1}{\hbar }\int _{0}^{t}{\rm d}s\;\int _{0}^{s}{\rm d}u\;\mathop{\sum }\limits_{\alpha =1}^{2}{{K}_{\alpha }}(u-s){{q}_{\alpha }}(s){{q}_{\alpha }}(u) \right\}. \\ \end{array}\end{eqnarray} \tag{ 9 }$$

with ${{Q}_{\alpha }}=\frac{1}{2}({{q}_{\alpha +}}+{{q}_{\alpha -}})$, ${{q}_{\alpha }}={{q}_{\alpha +}}-{{q}_{\alpha -}}$ and N(t) is a normalisation factor that can be determined from the conservation of the normalisation of the density matrix, ${\rm tr}{{\hat{\rho }}_{{\rm S}}}=1$,

$$\begin{eqnarray}&&N(t)={{\pi }^{2}}{{\hbar }^{2}}/\left| {{A}_{16}}(t){{A}_{38}}(t)-{{A}_{18}}(t){{A}_{36}}(t) \right|.\end{eqnarray} \tag{ 10 }$$

The matrix elements A_ij are defined in equation (A.49). The new coordinates ${{Q}_{\alpha }}$ and ${{q}_{\alpha }}$ satisfy the following equation of motion

$$\begin{eqnarray}\begin{array}{rcl} {{\ddot{Q}}_{1,2}}(s)+\omega _{1,2}^{2}{{Q}_{1,2}}(s)+\frac{c(s)}{{{m}_{1,2}}}{{Q}_{2,1}}(s)+\frac{{\rm d}}{{\rm d}s}\int _{0}^{s}{\rm d}u\;{{\gamma }_{1,2}}(s-u){{Q}_{1,2}}(u) & = & 0, \\ {{\ddot{q}}_{1,2}}(s)+\omega _{1,2}^{2}{{q}_{1,2}}(s)+\frac{c(s)}{{{m}_{1,2}}}{{q}_{2,1}}(s)-\frac{{\rm d}}{{\rm d}s}\int _{s}^{t}{\rm d}u\;{{\gamma }_{1,2}}(u-s){{q}_{1,2}}(u) & = & 0, \\ \end{array}\end{eqnarray} \tag{ 11 }$$

with the boundary conditions for ${{Q}_{\alpha }}(0)=Q_{\alpha }^{\prime }$, ${{Q}_{\alpha }}(t)=Q_{\alpha }^{\prime\prime }$, ${{q}_{\alpha }}(0)=q_{\alpha }^{\prime }$, ${{q}_{\alpha }}(t)=q_{\alpha }^{\prime\prime }$.

In the original frequency-and-mass degenerate formulation of the driving-assisted high-temperature entanglement scenario[4], the spectral densities ${{J}_{1}}(\omega )={{J}_{2}}(\omega )=J(\omega )$ were taken in the Ohmic form $J(\omega )=m\gamma \omega$ with γ denoting the standard coupling to the environment constant. This leads to local-in-time dissipation in equations (11) provided by the fact that ${{\gamma }_{1,2}}(s)=2\gamma \delta (s)$, and induces larger decay rates for the loss of entanglement (see below). To overcome this and to have a more general and complete characterisation of the dynamics that allow for a closer description of experimental realisation [6, 7], the influence of the bath on the system is characterised here by the spectral density

$$\begin{eqnarray}&&{{J}_{\alpha }}(\omega )={{m}_{\alpha }}{{\gamma }_{\alpha }}\omega \frac{\Omega _{\alpha }^{2}}{{{\omega }^{2}}+\Omega _{\alpha }^{2}},\end{eqnarray} \tag{ 12 }$$

where ${{\Omega }_{\alpha }}$ denotes a finite cut-off frequency. By replacing the last expression in equation (7), this spectral density generates ${{\gamma }_{\alpha }}(s)={{\gamma }_{\alpha }}{{\Omega }_{\alpha }}{\rm exp} \left( -{{\Omega }_{\alpha }}\left| s \right| \right),$ which leads to memory effects in the dissipation for times $s\lt {{\tau }_{\alpha }}=\Omega _{\alpha }^{-1}$ in the equations of motion (11). For this spectral density ${{K}_{\alpha }}(s)={{m}_{\alpha }}{{\gamma }_{\alpha }}\Omega _{\alpha }^{2}{{(\hbar {{\beta }_{\alpha }})}^{-1}}\sum _{n=-\infty }^{\infty }\left[ {{\Omega }_{\alpha }}{\rm exp} \left( -{{\Omega }_{\alpha }}|s| \right)-|{{\nu }_{\alpha ,n}}|{\rm exp} \left( -|{{\nu }_{\alpha ,n}}||s| \right) \right]{{\left[ \Omega _{\alpha }^{2}-\nu _{\alpha ,n}^{2} \right]}^{-1}},$ with ${{\nu }_{\alpha ,n}}=2\pi n{{(\hbar {{\beta }_{\alpha }})}^{-1}}$ being Matsubara's frequencies[26, 27, 32]. Note that under well-defined conditions, the spectral density can be experimentally reconstructed via linear absorption spectroscopy [38] and there is no need to assume a specific one.

For the degenerate Ohmic situation in Ref.[4], the solution to equation (11) can be expressed in terms of the solutions of Mathieu's oscillator (see also [39]). The non-Ohmic character of the spectral density in equation (12) prevents the formulation of the solution of equation (11) in a similar fashion. However, their linear character allows for expressing the formal solution in the form

$$\begin{eqnarray}\begin{array}{rcl} {{Q}_{1}}(t,s) & = & {{U}_{1}}(t,s)Q_{1}^{\prime }+{{U}_{2}}(t,s)Q_{1}^{\prime\prime }+{{U}_{3}}(t,s)Q_{2}^{\prime }+{{U}_{4}}(t,s)Q_{2}^{\prime\prime }, \\ {{Q}_{2}}(t,s) & = & {{V}_{1}}(t,s)Q_{2}^{\prime }+{{V}_{2}}(t,s)Q_{2}^{\prime\prime }+{{V}_{3}}(t,s)Q_{1}^{\prime }+{{V}_{4}}(t,s)Q_{1}^{\prime\prime }, \\ {{q}_{1}}(t,s) & = & {{u}_{1}}(t,s)q_{1}^{\prime }+{{u}_{2}}(t,s)q_{1}^{\prime\prime }+{{u}_{3}}(t,s)q_{2}^{\prime }+{{u}_{4}}(t,s)q_{2}^{\prime\prime }, \\ {{q}_{2}}(t,s) & = & {{v}_{1}}(t,s)q_{2}^{\prime }+{{v}_{2}}(t,s)q_{2}^{\prime\prime }+{{v}_{3}}(t,s)q_{1}^{\prime }+{{v}_{4}}(t,s)q_{1}^{\prime\prime }, \\ \end{array}\end{eqnarray} \tag{ 13 }$$

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 equation (11). Because of the time-reversed character of the limits in the integral contribution to the equation of motion of the coordinate ${{\hat{q}}_{\alpha }}$ in equation (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 equation (9).

Before proceeding, a note is in order. Strict non-Markovian dynamics are characterised by a flat (ω-independent) power noise $S(\omega )=\frac{1}{2\;m}J(\omega ){\rm coth} (\frac{1}{2}\hbar \omega \beta )$. Although the analytic formulae derived here are valid at any temperature regime, interest in this article is in the high temperature regime. Thus, for the temperature regime relevant to this article, ${\rm coth} (\frac{1}{2}\hbar \omega \beta )\sim 2{{k}_{{\rm B}}}T{{\hbar }^{-1}}{{\omega }^{-1}}$. Hence, for purely Ohmic dissipation $J(\omega )\sim \omega$, the power noise becomes ω-independent and the dynamics are Markovian. Thus, the case of purely Ohmic dissipation is referred to as Markovian dynamics whereas the case of Ohmic dissipation with a finite cut-off frequency is understood to account for non-Markovian dynamics. Further details on the power noise and its relation to the non-Markovian character of the dynamics can be found in [8, 36].

The logarithmic negativity is defined as ${{E}_{{\rm N}}}=-\frac{1}{2}\sum _{i=1}^{4}{{{\rm log} }_{2}}[{\rm min} (1,2|{{l}_{i}}|)],$ where l_i's are the symplectic eigenvalues of the covariance matrix. They are the normal eigenvalues of the matrix $-{\rm i}\Sigma \sigma$, with $\Sigma =\left( \begin{array}{ccccccccccccccc} 0 & {{\sf I}_{2}} \\ -{{\sf I}_{2}} & 0 \\ \end{array} \right)$ the symplectic matrix and ${{\sf I}_{2}}$ the identity matrix of dimension 2. For separable states, ${{\hat{\rho }}_{{\rm S}}}={{\sum }_{i}}{{p}_{i}}\hat{\rho }_{1}^{(i)}\;\otimes \;\hat{\rho }_{2}^{(i)}$, the logarithmic negativity of the system is zero and each oscillator can be described independently. For continuous variable systems, entanglement has as an upper limit the maximally entangled EPR wave function with ${{E}_{{\rm N}}}\to \infty .$ Hence, the amount of entanglement measured by the logarithmic negativity is unbounded from above.

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 with that oscillator. This is obtained by tracing out the reduced density matrix over the coordinates of the other oscillator; i.e., ${{\hat{\rho }}_{{{{\rm S}}_{1,2}}}}={\rm t}{{{\rm r}}_{{{{\rm S}}_{2,1}}}}{{\hat{\rho }}_{{\rm S}}}.$ For instance, the matrix elements of the reduced density operator associated with the first oscillator are given by

$$\begin{eqnarray}&&\left\langle Q_{1}^{\prime\prime } \right|{{\hat{\rho }}_{{{{\rm S}}_{1}}}}(t)\left| q_{1}^{\prime\prime } \right\rangle =\int _{-\infty }^{\infty }{\rm d}Q_{2}^{\prime\prime }\;\left\langle Q_{1}^{\prime\prime },Q_{2}^{\prime\prime } \right|{{\hat{\rho }}_{{\rm S}}}(t)\left| q_{1}^{\prime\prime },q_{2}^{\prime\prime }=0 \right\rangle ,\end{eqnarray} \tag{ 14 }$$

while for the second oscillator $\left\langle Q_{2}^{\prime\prime } \right|{{\hat{\rho }}_{{{{\rm S}}_{2}}}}(t)\left| q_{2}^{\prime\prime } \right\rangle =\int _{-\infty }^{\infty }{\rm d}Q_{1}^{\prime\prime }\;\left\langle Q_{1}^{\prime\prime },Q_{2}^{\prime\prime } \right|{{\hat{\rho }}_{{\rm S}}}(t)\left| q_{1}^{\prime\prime }=0,q_{2}^{\prime\prime } \right\rangle .$

To find the first and second moments of each oscillator, it is necessary to perform the following integrals:

$$\begin{eqnarray}\begin{array}{rcl} \left\langle {{{\hat{q}}}_{\alpha }} \right\rangle (t) & = & \int _{-\infty }^{\infty }{\rm d}Q_{\alpha }^{\prime\prime }\;Q_{\alpha }^{\prime\prime }\;\left\langle Q_{\alpha }^{\prime\prime } \right|{{{\hat{\rho }}}_{{{{\rm S}}_{\alpha }}}}(t)\left| q_{\alpha }^{\prime\prime }=0 \right\rangle , \\ \left\langle {{{\hat{p}}}_{\alpha }} \right\rangle (t) & = & -{\rm i}\hbar \int _{-\infty }^{\infty }{\rm d}Q_{\alpha }^{\prime\prime }\;\frac{{\rm d}}{{\rm d}q_{\alpha }^{\prime\prime }}{{\left. \left\langle Q_{\alpha }^{\prime\prime } \right|{{{\hat{\rho }}}_{{{{\rm S}}_{\alpha }}}}(t)\left| q_{\alpha }^{\prime\prime } \right\rangle \right|}_{q_{\alpha }^{\prime\prime }=0}}, \\ \left\langle \hat{q}_{\alpha }^{2} \right\rangle (t) & = & \int _{-\infty }^{\infty }{\rm d}Q_{\alpha }^{\prime\prime }\;Q_{\alpha }^{\prime\prime 2}\left\langle Q_{\alpha }^{\prime\prime } \right|{{{\hat{\rho }}}_{{{{\rm S}}_{\alpha }}}}(t)\left| q_{\alpha }^{\prime\prime }=0 \right\rangle , \\ \left\langle \hat{p}_{\alpha }^{2} \right\rangle (t) & = & -{{\hbar }^{2}}\int _{-\infty }^{\infty }{\rm d}Q_{\alpha }^{\prime\prime }\;\frac{{{{\rm d}}^{2}}}{{\rm d}q_{\alpha }^{{{^{\prime\prime} }^{2}}}}{{\left. \left\langle Q_{\alpha }^{\prime\prime } \right|{{{\hat{\rho }}}_{{{{\rm S}}_{\alpha }}}}(t)\left| q_{\alpha }^{\prime\prime } \right\rangle \right|}_{q_{\alpha }^{\prime\prime }=0}}, \\ \frac{1}{2}\left\langle {{{\hat{q}}}_{\alpha }}{{{\hat{p}}}_{\alpha }}+{{{\hat{p}}}_{\alpha }}{{{\hat{q}}}_{\alpha }} \right\rangle (t) & = & -{\rm i}\hbar \int _{-\infty }^{\infty }{\rm d}Q_{\alpha }^{\prime\prime }\;Q_{\alpha }^{\prime\prime }\frac{{\rm d}}{{\rm d}q_{\alpha }^{\prime\prime }}{{\left. \left\langle Q_{\alpha }^{\prime\prime } \right|{{{\hat{\rho }}}_{{{{\rm S}}_{\alpha }}}}(t)\left| q_{\alpha }^{\prime\prime } \right\rangle \right|}_{q_{\alpha }^{\prime\prime }=0}}. \\ \end{array}\end{eqnarray} \tag{ 15 }$$

To find the covariances between the position and/or momentum operators between the oscillators, it is necessary to integrate out

$$\begin{eqnarray}\begin{array}{rcl} \left\langle {{{\hat{q}}}_{1}}{{{\hat{q}}}_{2}} \right\rangle (t) & = & \int _{-\infty }^{\infty }{\rm d}Q_{1}^{\prime\prime }\;{\rm d}Q_{2}^{\prime\prime }\;Q_{1}^{\prime\prime }Q_{2}^{\prime\prime }\left\langle Q_{1}^{\prime\prime },Q_{2}^{\prime\prime } \right|{{{\hat{\rho }}}_{{\rm S}}}(t)\left| q_{1}^{\prime\prime }=0,q_{2}^{\prime\prime }=0 \right\rangle , \\ \left\langle {{{\hat{q}}}_{1}}{{{\hat{p}}}_{2}} \right\rangle (t) & = & -{\rm i}\hbar \int _{-\infty }^{\infty }{\rm d}Q_{1}^{\prime\prime }\;{\rm d}Q_{2}^{\prime\prime }\;Q_{1}^{\prime\prime }\frac{{\rm d}}{{\rm d}q_{2}^{\prime\prime }}{{\left. \left\langle Q_{1}^{\prime\prime },Q_{2}^{\prime\prime } \right|{{{\hat{\rho }}}_{{\rm S}}}(t)\left| q_{1}^{\prime\prime }=0,q_{2}^{\prime\prime } \right\rangle \right|}_{q_{2}^{\prime\prime }=0}}, \\ \left\langle {{{\hat{q}}}_{2}}{{{\hat{p}}}_{1}} \right\rangle (t) & = & -{\rm i}\hbar \int _{-\infty }^{\infty }{\rm d}Q_{1}^{\prime\prime }\;{\rm d}Q_{2}^{\prime\prime }\;Q_{2}^{\prime\prime }\frac{{\rm d}}{{\rm d}q_{1}^{\prime\prime }}{{\left. \left\langle Q_{1}^{\prime\prime },Q_{2}^{\prime\prime } \right|{{{\hat{\rho }}}_{{\rm S}}}(t)\left| q_{1}^{\prime\prime },q_{2}^{\prime\prime }=0 \right\rangle \right|}_{q_{1}^{\prime\prime }=0}}, \\ \left\langle {{{\hat{p}}}_{1}}{{{\hat{p}}}_{2}} \right\rangle (t) & = & -{{\hbar }^{2}}\int _{-\infty }^{\infty }{\rm d}Q_{1}^{\prime\prime }\;{\rm d}Q_{2}^{\prime\prime }\;\frac{{{{\rm d}}^{2}}}{{\rm d}q_{1}^{\prime\prime }{\rm d}q_{2}^{\prime\prime }}{{\left. \left\langle Q_{1}^{\prime\prime },Q_{2}^{\prime\prime } \right|{{{\hat{\rho }}}_{{\rm S}}}(t)\left| q_{1}^{\prime\prime },q_{2}^{\prime\prime } \right\rangle \right|}_{q_{1}^{\prime\prime }=0,q_{2}^{\prime\prime }=0}}. \\ \end{array}\end{eqnarray} \tag{ 16 }$$

For the case when no initial correlations exist between the oscillators in the system of interest, namely, ${{\hat{\rho }}_{{\rm S}}}(0)={{\hat{\rho }}_{{{{\rm S}}_{1}}}}(0)\;\otimes \;{{\hat{\rho }}_{{{{\rm 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 utilising the expressions (15) and (16) in the context, e.g., of optimal control of sideband cooling of nano-mechanical resonators [41].

To obtain specific results, both oscillators are assumed to be in a general Gaussian state, and therefore the initial density matrix for each oscillator can be expressed in terms of the coordinates ${{Q}_{\alpha }}$ and ${{q}_{\alpha }}$ as

$$\begin{eqnarray}\begin{array}{rcl} \langle Q_{\alpha }^{\prime }|{{{\hat{\rho }}}_{\alpha }}(0)|q_{\alpha }^{\prime }\rangle & = & \frac{1}{\sqrt{2\pi \sigma _{q}^{(\alpha )}}}{\rm exp} \left\{ -\frac{1}{2\sigma _{qq}^{(\alpha )}}{{\left( Q_{\alpha }^{\prime }-\left\langle {{q}_{\alpha }} \right\rangle \right)}^{2}}-\frac{1}{2{{\hbar }^{2}}}\left( \sigma _{pp}^{(\alpha )}-\frac{\sigma _{qp}^{2(\alpha )}}{\sigma _{qq}^{(\alpha )}} \right)q_{\alpha }^{\prime } \right. \\ {} & {} & +\left. \frac{{\rm i}}{\hbar }\left[ \left\langle {{p}_{\alpha }} \right\rangle +\frac{\sigma _{qp}^{(\alpha )}}{\sigma _{qq}^{(\alpha )}}\left( Q_{\alpha }^{\prime }-\left\langle {{q}_{\alpha }} \right\rangle \right) \right]q_{\alpha }^{\prime } \right\}, \\ \end{array}\end{eqnarray} \tag{ 17 }$$

where $\left\langle {{q}_{\alpha }} \right\rangle$, $\left\langle {{p}_{\alpha }} \right\rangle$, $\sigma _{qq}^{(\alpha )}$, $\sigma _{pp}^{(\alpha )}$ and $\sigma _{qp}^{(\alpha )}$ are the first moments of position and momentum, and the variance of the position, momentum and position-momentum, respectively, for the $\alpha {\rm th}$ oscillator.

Figure 2 shows the time evolution of the logarithmic negativity for a variety of coupling temperatures. Dark thick curves correspond to the non-Markovian case [$\Omega =\omega$ in equation (12)] and light thin curves correspond to the Markovian case [$\Omega =20\omega$ in equation (12)]. From there, it is clear that for the particular functional form of $J(\omega )$ in equation (12), non-Markovian dynamics are able to support the same amount of steady state entanglement at twice the temperature of the corresponding Markovian case.

An additional dynamic feature present in figures 1and 2 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 under non-Markovian dynamics (see section 4.4), it is then natural to expect that the driving force needs to perform less work to squeeze the normal modes of the oscillators when the dynamics are non-Markovian. Hence, if the power of the driving field is kept constant for the Markovian and the non-Markovian situations, entanglement will be generated faster in the non-Markovian case.

Figure 2 depicts the time dynamics of entanglement for a variety of coupling rates γ at fixed temperature. In analogy to the case discussed in figure 1, non-Markovian dynamics are able to support the same amount of steady state entanglement at twice the coupling rate γ of the corresponding Markovian case. Because of the non-Markovian character of the dynamics, the simulations over several periods of the driving force are very expensive in computational terms; it is therefore convenient to set the amplitude strength c₁ to a rather large value (${{c}_{1}}=0.2\;m{{\omega }^{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₁ (see section 4.4).

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 its initial state [4]. Despite the non-Markovian character of the dynamics and the associated dependence on the history of the system evolution, this feature remains present here. Figure 3 shows the time evolution of the logarithmic negativity for a variety of initial states given by equation (17). There, it is clear that all of them reach the same amount of stationary entanglement.

Zoom In Zoom Out Reset image size

If the characteristic frequency of a given system is denoted by ω, the survival of quantum features such as entanglement can be predicted for parameters that satisfy the condition $S(\omega )\lt \omega \mu$, with $S(\omega )=\frac{1}{2\;m}J(\omega ){\rm coth} (\frac{1}{2}\hbar \omega \beta )$ being the power noise and μ the pumping rate. Specifically, at high temperatures of $\hbar \omega \beta \ll 1$, entanglement survives in the steady state if

$$\begin{eqnarray}&&{{k}_{{\rm B}}}TJ(\omega )/m\lt \hbar {{\omega }^{2}}\mu .\end{eqnarray} \tag{ 19 }$$

Because of the various time/energy scales involved in this non-trivial out-of-equilibrium situation, the impact of non-Markovian dynamics can be 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 equation (12) is discussed below and three non-Markovianly scaled parameters, ${{T}_{{\rm nM}}}$, ${{\gamma }_{{\rm nM}}}$ and ${{\mu }_{{\rm nM}}}$, are introduced.

After plugging the spectral density (12) in the non-Markovian quantum limit given by equation (19), a new effective temperature can be defined

$$\begin{eqnarray}&&{{T}_{{\rm nM}}}=T\left( 1-\frac{1}{1+{{\Omega }^{2}}/{{\omega }^{2}}} \right),\end{eqnarray} \tag{ 20 }$$

such that the quantum limit can be cast in the form found in [4]; namely, ${{k}_{{\rm B}}}{{T}_{{\rm nM}}}/\hbar \omega \leqslant \mu /\gamma$, but with ${{T}_{{\rm nM}}}$ instead of T. It is also possible to define an effective coupling constant

$$\begin{eqnarray}&&{{\gamma }_{{\rm nM}}}=\gamma \left( 1-\frac{1}{1+{{\Omega }^{2}}/{{\omega }^{2}}} \right).\end{eqnarray} \tag{ 21 }$$

For $\Omega \sim \omega$, ${{T}_{{\rm nM}}}\sim \frac{1}{2}T$ and ${{\gamma }_{{\rm nM}}}\sim \frac{1}{2}\gamma$. This scaling of the temperature or the coupling constant explains the results depicted in figures 1 and 2, respectively. Alternatively, the non-Markovian scaling factor $1-\frac{1}{1+{{\Omega }^{2}}/{{\omega }^{2}}}$ can be assigned to a third energy scale. We define

$$\begin{eqnarray}&&{{\mu }_{{\rm nM}}}=\mu {{\left( 1-\frac{1}{1+{{\Omega }^{2}}/{{\omega }^{2}}} \right)}^{-1}},\end{eqnarray} \tag{ 22 }$$

so that ${{\mu }_{{\rm nM}}}\geqslant \mu$ provided by the fact that $1-\frac{1}{1+{{\Omega }^{2}}/{{\omega }^{2}}}\leqslant 1$. To be more concrete, note that the particular situation analysed in [4] corresponds to the case $J(\omega )/m=\gamma \omega$ and μ corresponds to the imaginary part of the associated Mathieu coefficient. Assuming that to leading order in the coupling c₁ [see equation (18)] the imaginary part of the Mathieu coefficient can still be expressed by $\mu \sim {{c}_{1}}/4\omega$, 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.

Although the non-Markovianly scaled parameters ${{T}_{{\rm nM}}}$, ${{\gamma }_{{\rm nM}}}$ and ${{\mu }_{{\rm nM}}}$ are particular to the spectral density in equation (12), based on the extensively studied features of non-Markovian dynamics [8, 10–12, 14–16, 19, 21, 22], there is no apparent physical reason not to expect the same scaling scenario in general.

From figure 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 temperatures is assisted here by the interaction between oscillators, the following should be considered: (i) the possible role of heat transfer in preserving/destroying quantum correlations between oscillators, and (ii) the classical/quantum nature of possible correlations established by heat fluxes at high/low temperatures among the oscillators. For concreteness of the present work, these concerns are analysed in a separated contribution and here our interest is restricted to the amount of entanglement that can be reached for different temperature ratios. Specifically, figure 4 depicts the time dynamics of entanglement for a variety of temperature ratios ${{T}_{2}}/{{T}_{1}}$ at fixed T₁. The main features in this figure are (i) the strong dependence of the time at which entanglement is generated on the temperature ratio and (ii) 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.

Zoom In Zoom Out Reset image size

The effective temperature at which each oscillator thermalises is a function of the power noise of the bath, and therefore of the coupling constant [8]. Thus, another interesting situation from a thermodynamic point of view is the case when the coupling constants are different. Figure 5 depicts the time dynamics of entanglement for a variety of coupling constant ratios ${{\gamma }_{2}}/{{\gamma }_{1}}$ at fixed ${{\gamma }_{1}}$. Up to the weak dependence of the time at which entanglement is generated on the coupling ratios, the behaviour is essentially the same as in figure 4. Below, the case of non-degenerate oscillators is analysed and the effect of non-Markovian dynamics in the resonance condition ${{\omega }_{{\rm d}}}={{\omega }_{1}}+{{\omega }_{2}}$ is discussed.

Zoom In Zoom Out Reset image size

For small values of c₁ and in the rotating wave approximation, the maximum rate of generation of squeezing is obtained for ${{\omega }_{{\rm d}}}={{\omega }_{1}}+{{\omega }_{2}}$. For the degenerate case, it reads ${{\omega }_{{\rm d}}}=2{{\omega }_{1}}$. Below, for this optimal condition, the dynamics of entanglement are analysed for a variety of frequency ratios ${{\omega }_{2}}/{{\omega }_{1}}$ for fixed ${{\omega }_{1}}$. To isolate the effect of non-Markovian dynamics, parameters are chosen so that no entanglement is found in the Markovian degenerate case with ${{\omega }_{{\rm d}}}=2{{\omega }_{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 ($\sim 7\%$) and with 'non-resonant' driving. This feature adds to the known robustness of the squeezing-entanglement generation against small detuning from the resonance frequency.

Zoom In Zoom Out Reset image size

In the experimental work on the generation of two-mode squeezing reported in [6], ${{\omega }_{2}}\sim 0.939{{\omega }_{1}}$ with ${{\omega }_{{\rm d}}}={{\omega }_{1}}+{{\omega }_{2}}$ and $\gamma \sim {{10}^{-3}}$, so that for the non-Markovian scenario and pumping rate in figure 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₁. As 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 idea of the effectiveness of the driving field in the presence of lighter masses.

Zoom In Zoom Out Reset image size

To derive the influence functional, consider the following property of the influence functional:

If a system is interacting simultaneously with two uncoupled and independent environments A and B, with no initial correlations between them, then [23]

$$\begin{eqnarray}&&F={{F}_{A}}\cdot {{F}_{B}}.\end{eqnarray} \tag{ A.15 }$$

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,

$$\begin{eqnarray}&&\mathcal{F}[{{q}_{1+}},{{q}_{2+}},{{q}_{1-}},{{q}_{2-}}]=\mathcal{F}[{{q}_{1+}},{{q}_{1-}}]\mathcal{F}[{{q}_{2+}},{{q}_{2-}}],\end{eqnarray} \tag{ A.16 }$$

where

$$\begin{eqnarray}\begin{array}{rcl} \mathcal{F}[{{q}_{\alpha +}},{{q}_{\alpha -}}] & = & \int _{-\infty }^{\infty }{\rm d}{\boldsymbol{q}} _{\alpha +}^{\prime\prime }\;{\rm d}{\boldsymbol{q}} _{\alpha +}^{\prime }\;{\rm d}{\boldsymbol{q}} _{\alpha -}^{\prime }\;{{\rho }_{{{{\rm B}}_{\alpha }}}}({\boldsymbol{q}} _{\alpha +}^{\prime },{\boldsymbol{q}} _{\alpha -}^{\prime },0) \\ {} & {} & \times \int \mathcal{D}{{{\boldsymbol{q}} }_{\alpha +}}\;\mathcal{D}{{{\boldsymbol{q}} }_{\alpha -}}\;{\rm exp} \left\{ \frac{{\rm i}}{\hbar }\left( {{S}_{{{{\rm B}}_{\alpha }}}}[{{{\boldsymbol{q}} }_{\alpha +}}]+{{S}_{{{{\rm I}}_{\alpha }}}}[{{{\boldsymbol{q}} }_{\alpha +}},{{{\boldsymbol{q}} }_{\alpha +}}] \right. \right. \\ {} & {} & -\left. \left. {{S}_{{{{\rm B}}_{\alpha }}}}[{{{\boldsymbol{q}} }_{\alpha -}}]-{{S}_{{{{\rm I}}_{\alpha }}}}[{{q}_{\alpha -}},{{{\boldsymbol{q}} }_{\alpha -}}] \right) \right\},\quad \alpha =1,2. \\ \end{array}\end{eqnarray} \tag{ A.17 }$$

Having in mind that the oscillators encompassing each thermal bath are independent of them, the property of the influence functional stated above can be used once more so that each influence functional in equation (A.16) is given by

$$\begin{eqnarray}&&\mathcal{F}[{{q}_{\alpha +}},{{q}_{\alpha -}}]=\mathop{\prod }\limits_{k=1}^{N}{{\mathcal{F}}_{k}}[{{q}_{\alpha +}},{{q}_{\alpha -}}],\end{eqnarray} \tag{ A.18 }$$

where each ${{\mathcal{F}}_{k}}[{{q}_{\alpha +}},{{q}_{\alpha -}}]$ describes the influence of each oscillator in one thermal bath on the oscillator in the system of interest. Each of these influence functionals reads

$$\begin{eqnarray}\begin{array}{rcl} {{\mathcal{F}}_{k}}[q_{\alpha +}^{\prime\prime },{{q}_{\alpha -}}] & = & \int _{-\infty }^{\infty }{\rm d}{{q}_{\alpha ,k+}}\;{\rm d}q_{\alpha ,k+}^{\prime }\;{\rm d}q_{\alpha ,k-}^{\prime }\;\langle q_{\alpha ,k+}^{\prime }|\rho _{{{{\rm B}}_{\alpha }}}^{(k)}(0)|q_{\alpha ,k-}^{\prime }\rangle \\ {} & {} & \times \int \mathcal{D}{{q}_{\alpha ,k+}}\;{\rm exp} \left\{ \frac{{\rm i}}{\hbar }\left( S_{{{{\rm B}}_{\alpha }}}^{(k)}[{{q}_{\alpha ,k+}}]+S_{{{{\rm I}}_{\alpha }}}^{(k)}[{{q}_{\alpha +}},{{q}_{\alpha ,k+}}] \right) \right\} \\ {} & {} & \times \int \mathcal{D}{{q}_{\alpha ,k-}}\;{\rm exp} \left\{ -\frac{{\rm i}}{\hbar }\left( S_{{{{\rm B}}_{\alpha }}}^{(k)}[{{q}_{\alpha ,k-}}]+S_{{{{\rm I}}_{\alpha }}}^{(k)}[{{q}_{\alpha -}},{{q}_{\alpha ,k-}}] \right) \right\}. \\ \end{array}\end{eqnarray} \tag{ A.19 }$$

$\langle q_{\alpha ,k+}^{\prime }|\rho _{{{{\rm B}}_{\alpha }}}^{(k)}(0)|q_{\alpha ,k-}^{\prime }\rangle$ denotes the initial density matrix for the $k{\rm th}$ oscillator in the thermal bath. For the case of a thermal bath at thermal equilibrium at a temperature ${{T}_{\alpha }}$, it reads (see, e.g., [32])

$$\begin{eqnarray}&&\begin{array}{l} \langle q_{\alpha ,k+}^{\prime }|\rho _{{{{\rm B}}_{\alpha }}}^{(k)}(0)|q_{\alpha ,k-}^{\prime }\rangle \\ \quad =2{\rm sinh} \left( \frac{\hbar {{\beta }_{\alpha }}{{\omega }_{\alpha ,k}}}{2} \right)\sqrt{\frac{{{m}_{\alpha ,k}}{{\omega }_{\alpha ,k}}}{2\pi \hbar {\rm sinh} (\hbar {{\beta }_{\alpha }}{{\omega }_{\alpha ,k}})}} \\ \qquad \times {\rm exp} \left\{ -\frac{{{m}_{\alpha ,k}}{{\omega }_{\alpha ,k}}}{2\hbar {\rm sinh} (\hbar {{\beta }_{\alpha }}{{\omega }_{\alpha ,k}})}\left[ \left( q_{\alpha ,k+}^{\prime 2}+q_{\alpha ,k-}^{\prime 2} \right){\rm cosh} (\hbar {{\beta }_{\alpha }}{{\omega }_{\alpha ,k}})-2q_{\alpha ,k+}^{\prime }q_{\alpha ,k-}^{\prime } \right] \right\}. \\ \end{array}\end{eqnarray} \tag{ A.20 }$$

The next step is to evaluate the path integrals in equation (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_{{{{\rm B}}_{\alpha }}}^{(k)}$ and $S_{{{{\rm I}}_{\alpha }}}^{(k)}$ along the classical paths ${{q}_{\alpha ,k}}(s)$ of each oscillator in the thermal baths. These classical paths are the solution to the equation of motion

$$\begin{eqnarray}&&{{\ddot{q}}_{\alpha ,k}}(s)+\omega _{\alpha ,k}^{2}{{q}_{\alpha ,k}}(s)={{q}_{\alpha }}(s)\frac{{{c}_{\alpha ,k}}}{{{m}_{\alpha ,k}}},\end{eqnarray} \tag{ A.21 }$$

that can be obtained from the Lagrangians in (A.7) and (A.8). The trick for solving this differential equation lies in treating the system coordinate ${{q}_{\alpha }}$ 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., [32]), except for the term

$$\begin{eqnarray*}&&q_{\alpha }^{2}\frac{c_{\alpha ,k}^{2}}{2{{m}_{\alpha ,k}}\omega _{\alpha ,k}^{2}}.\end{eqnarray*}$$

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 equation (A.4), taking ${{{\boldsymbol{q}} }_{\alpha -}}(t)={\boldsymbol{q}} _{\alpha +}^{\prime\prime }$ for the tracing operation and using the propagation kernel for a driven harmonic oscillator, the path integrals in (A.19) are readily given by

$$\begin{eqnarray}&&\begin{array}{l} \int \mathcal{D}{{q}_{\alpha ,k\pm }}\;{\rm exp} \left\{ \pm \frac{{\rm i}}{\hbar }\left( S_{{{{\rm B}}_{\alpha }}}^{(k)}[{{q}_{\alpha ,k\pm }}]+S_{{{{\rm I}}_{\alpha }}}^{(k)}[{{q}_{\alpha \pm }},{{q}_{\alpha ,k\pm }}] \right) \right\} \\ \quad =\sqrt{\frac{{{m}_{\alpha ,k}}{{\omega }_{\alpha ,k}}}{2\pi (\pm {\rm i})\hbar {\rm sin} ({{\omega }_{\alpha ,k}}t)}}{\rm exp} \left( \pm \frac{{\rm i}}{\hbar }\left\{ \frac{{{m}_{\alpha ,k}}{{\omega }_{\alpha ,k}}}{2{\rm sin} ({{\omega }_{\alpha ,k}}t)}\left[ \left( q_{\alpha ,k+}^{\prime\prime 2}+q_{\alpha ,k\pm }^{\prime 2} \right){\rm cos} ({{\omega }_{\alpha ,k}}t) \right. \right. \right. \\ \qquad -\left. 2q_{\alpha ,k+}^{\prime\prime }q_{\alpha ,k\pm }^{\prime } \right]+\frac{{{c}_{\alpha ,k}}q_{\alpha ,k+}^{\prime\prime }}{{\rm sin} ({{\omega }_{\alpha ,k}}t)}\int _{0}^{t}{\rm d}s\;{\rm sin} ({{\omega }_{\alpha ,k}}s){{q}_{\alpha \pm }}(s) \\ \qquad +\frac{{{c}_{\alpha ,k}}q_{\alpha ,k\pm }^{\prime }}{{\rm sin} ({{\omega }_{\alpha ,k}}t)}\int _{0}^{t}{\rm d}s\;{\rm sin} [{{\omega }_{\alpha ,k}}(t-s)]{{q}_{\alpha \pm }}(s)-\frac{c_{\alpha ,k}^{2}}{2{{m}_{\alpha ,k}}\omega _{\alpha ,k}^{2}}\int _{0}^{t}{\rm d}s\;q_{\alpha \pm }^{2}(s) \\ \qquad -\left. \left. \frac{c_{\alpha ,k}^{2}}{{{m}_{\alpha ,k}}{{\omega }_{\alpha ,k}}{\rm sin} ({{\omega }_{\alpha ,k}}t)}\int _{0}^{t}{\rm d}s\;\int _{0}^{s}{\rm d}u\;{\rm sin} ({{\omega }_{\alpha ,k}}u){\rm sin} [{{\omega }_{\alpha ,k}}(t-s)]{{q}_{\alpha \pm }}(s){{q}_{\alpha \pm }}(u) \right\} \right). \\ \end{array}\end{eqnarray} \tag{ A.22 }$$

Using this expression and equation (A.20), the influence functional in equation (A.19) takes the form

$$\begin{eqnarray}&&\begin{array}{l} \mathcal{F}[{{q}_{\alpha ,+}},{{q}_{\alpha ,-}}] \\ \quad =2{\rm sinh} \left( \frac{\hbar {{\beta }_{\alpha }}{{\omega }_{\alpha ,k}}}{2} \right)\sqrt{\frac{{{m}_{\alpha ,k}}{{\omega }_{\alpha ,k}}}{2\pi \hbar {\rm sinh} (\hbar {{\beta }_{\alpha }}{{\omega }_{\alpha ,k}})}}\sqrt{\frac{{{m}_{\alpha ,k}}{{\omega }_{\alpha ,k}}}{2\pi {\rm i}\hbar {\rm sin} ({{\omega }_{\alpha ,k}}t)}}\sqrt{\frac{{{m}_{\alpha ,k}}{{\omega }_{\alpha ,k}}{\rm i}}{2\pi \hbar {\rm sin} ({{\omega }_{\alpha ,k}}t)}} \\ \qquad \times \int _{-\infty }^{\infty }{\rm d}q_{\alpha ,k+}^{\prime\prime }\;{\rm d}q_{\alpha ,k+}^{\prime }\;{\rm d}q_{\alpha ,k-}^{\prime }\;{\rm exp} \left\{ -\frac{{{m}_{\alpha ,k}}{{\omega }_{\alpha ,k}}}{2\hbar {\rm sinh} (\hbar {{\beta }_{\alpha }}{{\omega }_{\alpha ,k}})}\left[ \left( q_{\alpha ,k+}^{\prime 2}+q_{\alpha ,k-}^{\prime 2} \right) \right. \right. \\ \qquad \times \left. {\rm cosh} (\hbar {{\beta }_{\alpha }}{{\omega }_{\alpha ,k}})-2q_{\alpha ,k+}^{\prime }q_{\alpha ,k-}^{\prime } \right]+\frac{{\rm i}}{\hbar }\left\{ \frac{{{m}_{\alpha ,k}}{{\omega }_{\alpha ,k}}}{2{\rm sin} ({{\omega }_{\alpha ,k}}t)}\left[ \left( q_{\alpha ,k+}^{\prime 2}-q_{\alpha ,k-}^{\prime 2} \right){\rm cos} ({{\omega }_{\alpha ,k}}t) \right. \right. \\ \qquad -\left. 2q_{\alpha ,k+}^{\prime\prime }(q_{\alpha ,k+}^{\prime }-q_{\alpha ,k-}^{\prime }) \right]+\frac{{{c}_{\alpha ,k}}{{q}_{\alpha ,k+^{\prime} }}}{{\rm sin} ({{\omega }_{\alpha ,k}}t)}\int _{0}^{t}{\rm d}s\;{\rm sin} [{{\omega }_{\alpha ,k}}(t-s)]{{q}_{\alpha +}}(s) \\ \qquad -\frac{{{c}_{\alpha ,k}}q_{\alpha ,k-}^{\prime }}{{\rm sin} ({{\omega }_{\alpha ,k}}t)}\int _{0}^{t}{\rm d}s\;{\rm sin} [{{\omega }_{\alpha ,k}}(t-s)]{{q}_{\alpha -}}(s)-\frac{c_{\alpha ,k}^{2}}{2{{m}_{\alpha ,k}}\omega _{\alpha ,k}^{2}}\int _{0}^{t}{\rm d}s\;\left[ q_{\alpha +}^{2}(s)-q_{\alpha -}^{2}(s) \right] \\ \qquad +\frac{{{c}_{\alpha ,k}}q_{\alpha ,k+}^{\prime\prime }}{{\rm sin} ({{\omega }_{\alpha ,k}}t)}\int _{0}^{t}{\rm d}s\;{\rm sin} ({{\omega }_{\alpha ,k}}s)\left[ {{q}_{\alpha +}}(s)-{{q}_{\alpha -}}(s) \right]-\frac{c_{\alpha ,k}^{2}}{{{m}_{\alpha ,k}}{{\omega }_{\alpha ,k}}{\rm sin} ({{\omega }_{\alpha ,k}}t)} \\ \qquad \times \left. \int _{0}^{t}{\rm d}s\;\int _{0}^{s}{\rm d}u\;{\rm sin} ({{\omega }_{\alpha ,k}}u){\rm sin} [{{\omega }_{\alpha ,k}}(t-s)]\left[ {{q}_{\alpha +}}(s){{q}_{\alpha +}}(u)-{{q}_{\alpha -}}(s){{q}_{\alpha -}}(u) \right] \right\}. \\ \end{array}\end{eqnarray} \tag{ A.23 }$$

After performing the integrations over $q_{\alpha ,k+}^{\prime\prime }$, $q_{\alpha ,k+}^{\prime }$, $q_{\alpha ,k-}^{\prime }$, the influence functional reads

$$\begin{eqnarray}&&\begin{array}{l} {{\mathcal{F}}_{k}}[{{q}_{\alpha +}},{{q}_{\alpha -}}] \\ \quad ={\rm exp} \left( -\frac{{\rm i}}{\hbar }\left\{ (q_{\alpha +}^{\prime }+q_{\alpha -}^{\prime })\int _{0}^{t}{\rm d}s\;\frac{c_{\alpha ,k}^{2}}{2{{m}_{\alpha ,k}}\omega _{\alpha ,k}^{2}}{\rm cos} ({{\omega }_{\alpha ,k}}s)\left[ {{q}_{\alpha +}}(s)-{{q}_{\alpha -}}(s) \right] \right. \right. \\ \qquad +\left. \left. \int _{0}^{t}{\rm d}s\;\int _{0}^{s}{\rm d}u\;\frac{c_{\alpha ,k}^{2}}{2{{m}_{\alpha ,k}}\omega _{\alpha ,k}^{2}}{\rm cos} [{{\omega }_{\alpha ,k}}(s-u)]\left[ {{{\dot{q}}}_{\alpha +}}(u)+{{{\dot{q}}}_{\alpha -}}(u) \right]\left[ {{q}_{\alpha +}}(s)-{{q}_{\alpha -}}(s) \right] \right\} \right) \\ \qquad \times {\rm exp} \left\{ -\frac{1}{\hbar }\frac{c_{\alpha ,k}^{2}}{2{{m}_{\alpha ,k}}{{\omega }_{\alpha ,k}}}\int _{0}^{t}{\rm d}s\;\int _{0}^{s}{\rm d}u\;{\rm coth} \left( \frac{\hbar {{\beta }_{\alpha }}{{\omega }_{\alpha ,k}}}{2} \right){\rm cos} [{{\omega }_{\alpha ,k}}(u-s)] \right. \\ \qquad \times \left[ {{q}_{\alpha +}}(s)-{{q}_{\alpha -}}(s) \right]\left[ {{q}_{\alpha +}}(u)-{{q}_{\alpha -}}(u) \right]\left. {} \right\}. \\ \end{array}\end{eqnarray} \tag{ A.24 }$$

The above expression describes the influence of the $k{\rm th}$ oscillator in the $\alpha {\rm th}$ thermal bath on the $\alpha {\rm 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,

$$\begin{eqnarray}&&\begin{array}{l} \mathcal{F}[{{q}_{\alpha +}},{{q}_{\alpha -}}] \\ \quad ={\rm exp} \left( -\frac{{\rm i}}{\hbar }\left\{ (q_{\alpha +}^{\prime }+q_{\alpha -}^{\prime })\int _{0}^{t}{\rm d}s\;\mathop{\sum }\limits_{k=1}^{N}\frac{c_{\alpha ,k}^{2}}{2{{m}_{\alpha ,k}}\omega _{\alpha ,k}^{2}}{\rm cos} ({{\omega }_{\alpha ,k}}s)\left[ {{q}_{\alpha +}}(s)-{{q}_{\alpha -}}(s) \right] \right. \right. \\ \qquad +\int _{0}^{t}{\rm d}s\;\int _{0}^{s}{\rm d}u\;\mathop{\sum }\limits_{k=1}^{N}\frac{c_{\alpha ,k}^{2}}{2{{m}_{\alpha ,k}}\omega _{\alpha ,k}^{2}}{\rm cos} [{{\omega }_{\alpha ,k}}(s-u)]\left[ {{{\dot{q}}}_{\alpha +}}(u)+{{{\dot{q}}}_{\alpha -}}(u) \right] \\ \qquad \times \left. \left. \left[ {{q}_{\alpha +}}(s)-{{q}_{\alpha -}}(s) \right] \right\} \right)\times {\rm exp} \left\{ -\frac{1}{\hbar }\int _{0}^{t}{\rm d}s\;\int _{0}^{s}{\rm d}u\;\mathop{\sum }\limits_{k=1}^{N}\frac{c_{\alpha ,k}^{2}}{2{{m}_{\alpha ,k}}{{\omega }_{\alpha ,k}}}{\rm coth} \left( \frac{\hbar {{\beta }_{\alpha }}{{\omega }_{\alpha ,k}}}{2} \right) \right. \\ \qquad \times \left. {\rm cos} [{{\omega }_{\alpha ,k}}(u-s)]\left[ {{q}_{\alpha +}}(s)-{{q}_{\alpha -}}(s) \right]\left[ {{q}_{\alpha +}}(u)-{{q}_{\alpha -}}(u) \right] \right\}. \\ \end{array}\end{eqnarray} \tag{ A.25 }$$

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}_{\alpha }}({{\omega }_{\alpha }})$, which is defined as

$$\begin{eqnarray}&&{{J}_{\alpha }}({{\omega }_{\alpha }})=\pi \mathop{\sum }\limits_{k=1}^{N}\frac{c_{\alpha ,k}^{2}}{2{{m}_{\alpha ,k}}{{\omega }_{\alpha ,k}}}\delta ({{\omega }_{\alpha }}-{{\omega }_{\alpha ,k}}).\end{eqnarray} \tag{ A.26 }$$

Therefore, for a discrete set of 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 ${{\omega }_{\alpha ,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 ${{\omega }_{\alpha }}$ [26]. This allows for defining the dissipation ${{\gamma }_{\alpha }}(s)$ and noise ${{K}_{\alpha }}(s)$ kernels in terms of the spectral density as

$$\begin{eqnarray}&&{{\gamma }_{\alpha }}(s)=\frac{2}{{{m}_{\alpha }}}\int _{0}^{\infty }\frac{{\rm d}{{\omega }_{\alpha }}\;}{\pi }\frac{{{J}_{\alpha }}({{\omega }_{\alpha }})}{{{\omega }_{\alpha }}}{\rm cos} ({{\omega }_{\alpha }}s),\end{eqnarray} \tag{ A.27 }$$

$$\begin{eqnarray}&&{{K}_{\alpha }}(s)=\int _{0}^{\infty }\frac{{\rm d}{{\omega }_{\alpha }}\;}{\pi }{{J}_{\alpha }}({{\omega }_{\alpha }}){\rm coth} \left( \frac{\hbar {{\beta }_{\alpha }}{{\omega }_{\alpha }}}{2} \right){\rm cos} ({{\omega }_{\alpha }}s).\end{eqnarray} \tag{ A.28 }$$

By replacing (A.26) in the last expressions,

$$\begin{eqnarray}&&{{\gamma }_{\alpha }}(s)=\frac{1}{{{m}_{\alpha }}}\mathop{\sum }\limits_{k=1}^{N}\frac{c_{\alpha ,k}^{2}}{2{{m}_{\alpha ,k}}\omega _{\alpha ,k}^{2}}{\rm cos} ({{\omega }_{\alpha ,k}}s),\end{eqnarray} \tag{ A.29 }$$

$$\begin{eqnarray}&&{{K}_{\alpha }}(s)=\mathop{\sum }\limits_{k=1}^{N}\frac{c_{\alpha ,k}^{2}}{2{{m}_{\alpha ,k}}{{\omega }_{\alpha ,k}}}{\rm coth} \left( \frac{\hbar {{\beta }_{\alpha }}{{\omega }_{\alpha ,k}}}{2} \right){\rm cos} ({{\omega }_{\alpha ,k}}s),\end{eqnarray} \tag{ A.30 }$$

which can be identified with the summations in the argument of the exponentials in equation (A.25). Having in mind these considerations, and using the expressions (A.26), (A.27) and (A.28), equation (A.25) can be expressed in terms of the dissipation and noise kernels as

$$\begin{eqnarray}&&\begin{array}{l} \mathcal{F}[{{q}_{\alpha +}},{{q}_{\alpha -}}] \\ \quad ={\rm exp} \left( -\frac{{\rm i}}{\hbar }\frac{{{m}_{\alpha }}}{2}\left\{ (q_{\alpha +}^{\prime }+q_{\alpha -}^{\prime })\int _{0}^{t}{\rm d}s\;{{\gamma }_{\alpha }}(s)\left[ {{q}_{\alpha +}}(s)-{{q}_{\alpha -}}(s) \right] \right. \right. \\ \qquad +\left. \left. \int _{0}^{t}{\rm d}s\;\int _{0}^{s}{\rm d}u\;{{\gamma }_{\alpha }}(s-u)\left[ {{{\dot{q}}}_{\alpha +}}(u)+{{{\dot{q}}}_{\alpha -}}(u) \right]\left[ {{q}_{\alpha +}}(s)-{{q}_{\alpha -}}(u) \right] \right\} \right) \\ \qquad \times {\rm exp} \left\{ -\frac{1}{\hbar }\int _{0}^{t}{\rm d}s\;\int _{0}^{s}{\rm d}u\;{{K}_{\alpha }}(u-s)\left[ {{q}_{\alpha +}}(s)-{{q}_{\alpha -}}(s) \right]\left[ {{q}_{\alpha +}}(u)-{{q}_{\alpha -}}(u) \right] \right\}. \\ \end{array}\end{eqnarray} \tag{ A.31 }$$

By replacing the above expression in (A.16), the total influence functional for the system of interest reads

$$\begin{eqnarray}&&\begin{array}{l} \mathcal{F}[{{q}_{1+}},{{q}_{2+}},{{q}_{1-}},{{q}_{2-}}] \\ \quad =\mathop{\prod }\limits_{\alpha =1}^{2}{\rm exp} \left( -\frac{{\rm i}}{\hbar }\frac{{{m}_{\alpha }}}{2}\left\{ (q_{\alpha +}^{\prime }+q_{\alpha -}^{\prime })\int _{0}^{t}{\rm d}s\;{{\gamma }_{\alpha }}(s)\left[ {{q}_{\alpha +}}(s)-{{q}_{\alpha -}}(s) \right] \right. \right. \\ \qquad +\left. \left. \int _{0}^{t}{\rm d}s\;\int _{0}^{s}{\rm d}u\;{{\gamma }_{\alpha }}(s-u)\left[ {{{\dot{q}}}_{\alpha +}}(u)+{{{\dot{q}}}_{\alpha -}}(u) \right]\left[ {{q}_{\alpha +}}(s)-{{q}_{\alpha -}}(u) \right] \right\} \right) \\ \qquad \times {\rm exp} \left\{ -\frac{1}{\hbar }\int _{0}^{t}{\rm d}s\;\int _{0}^{s}{\rm d}u\;{{K}_{\alpha }}(u-s)\left[ {{q}_{\alpha +}}(s)-{{q}_{\alpha -}}(s) \right]\left[ {{q}_{\alpha +}}(u)-{{q}_{\alpha -}}(u) \right] \right\}. \\ \end{array}\end{eqnarray} \tag{ A.32 }$$

To provide an explicit expression for the propagating function associated with the system of interest, consider the expression (A.13)

$$\begin{eqnarray}&&\begin{array}{l} J(q_{1+}^{\prime\prime },q_{2+}^{\prime\prime },q_{1-}^{\prime\prime },q_{2-}^{\prime\prime },t;q_{1+}^{\prime },q_{2+}^{\prime },q_{1-}^{\prime },q_{2-}^{\prime },0) \\ \quad =\int \mathcal{D}{{q}_{1+}}\;\mathcal{D}{{q}_{2+}}\;\mathcal{D}{{q}_{1-}}\;\mathcal{D}{{q}_{2-}}\;{\rm exp} \left\{ \frac{{\rm i}}{\hbar }\left( {{S}_{{\rm S}}}[{{q}_{1+}},{{q}_{2+}}]-{{S}_{{\rm S}}}[{{q}_{1-}},{{q}_{2-}}] \right) \right\} \\ \qquad \times \mathcal{F}[{{q}_{1+}},{{q}_{2+}},{{q}_{1-}},{{q}_{2-}}], \\ \end{array}\end{eqnarray} \tag{ A.33 }$$

with $\mathcal{F}[{{q}_{1+}},{{q}_{2+}},{{q}_{1-}},{{q}_{2-}}]$ given by the expression (A.32). Since the path integrals in (A.33) are quadratic, they can be performed exactly[37]. This leaves for the propagating function [4, 29, 30, 39] the following expression

$$\begin{eqnarray}&&\begin{array}{l} J(q_{1+}^{\prime\prime },q_{2+}^{\prime\prime },q_{1-}^{\prime\prime },q_{2-}^{\prime\prime },t;q_{1+}^{\prime },q_{2+}^{\prime },q_{1-}^{\prime },q_{2-}^{\prime },0) \\ \quad =\frac{1}{N(t)}{\rm exp} \left\{ \frac{{\rm i}}{\hbar }\left( {{S}_{{\rm S}}}[{{{\bar{q}}}_{1+}},{{{\bar{q}}}_{2+}}]-{{S}_{{\rm S}}}[{{{\bar{q}}}_{1-}},{{{\bar{q}}}_{2-}}] \right) \right\}\mathcal{F}[{{{\bar{q}}}_{1+}},{{{\bar{q}}}_{2+}},{{{\bar{q}}}_{1-}},{{{\bar{q}}}_{2-}}], \\ \end{array}\end{eqnarray} \tag{ A.34 }$$

where ${{\bar{q}}_{\alpha \pm }}$ denotes the classical paths of the system and N(t) is a normalisation factor so that ${\rm tr}{{\hat{\rho }}_{{\rm S}}}(t)=1$. Taking into account expressions (A.5) to (A.8) and (A.32), the propagating function takes the form

$$\begin{eqnarray}&&\begin{array}{l} J(q_{1+}^{\prime\prime },q_{2+}^{\prime\prime },q_{1-}^{\prime\prime },q_{2-}^{\prime\prime },t;q_{1+}^{\prime },q_{2+}^{\prime },q_{1-}^{\prime },q_{2-}^{\prime },0) \\ \quad =\frac{1}{N(t)}{\rm exp} \left\{ \frac{{\rm i}}{\hbar }\int _{0}^{t}{\rm d}s\;\left[ \mathop{\sum }\limits_{\alpha =1}^{2}\left( \frac{1}{2}{{m}_{\alpha }}\dot{{\bar{q}}}_{\alpha +}^{2}(s)-\frac{1}{2}{{m}_{\alpha }}\omega _{\alpha }^{2}\bar{q}_{\alpha +}^{2}(s) \right)-c(s){{{\bar{q}}}_{1+}}(s){{{\bar{q}}}_{2+}}(s) \right] \right\} \\ \qquad \times {\rm exp} \left\{ -\frac{{\rm i}}{\hbar }\int _{0}^{t}{\rm d}s\;\left[ \mathop{\sum }\limits_{\alpha =1}^{2}\left( \frac{1}{2}{{m}_{\alpha }}\dot{{\bar{q}}}_{\alpha -}^{2}(s)-\frac{1}{2}{{m}_{\alpha }}\omega _{\alpha }^{2}\bar{q}_{\alpha -}^{2}(s) \right)-c(s){{{\bar{q}}}_{1-}}(s){{{\bar{q}}}_{2-}}(s) \right] \right\} \\ \qquad \times \mathcal{F}[{{{\bar{q}}}_{1+}},{{{\bar{q}}}_{2+}},{{{\bar{q}}}_{1-}},{{{\bar{q}}}_{2-}}]. \\ \end{array}\end{eqnarray} \tag{ A.35 }$$

To simplify the last expression, consider a new set of coordinates[4, 29, 39]

$$\begin{eqnarray}&&{{Q}_{\alpha }}=\frac{1}{2}({{\bar{q}}_{\alpha +}}+{{\bar{q}}_{\alpha -}}),\quad {{q}_{\alpha }}={{\bar{q}}_{\alpha +}}-{{\bar{q}}_{\alpha -}},\end{eqnarray} \tag{ A.36 }$$

where the Jacobian of the transformation is equal to one. Hence, equation (A.35) takes the form

$$\begin{eqnarray}&&\begin{array}{l} J({\bf Q}^{{\prime} ^{\prime}} ,{\bf q}^{{\prime} ^{\prime}} ,t;{\bf Q}^{\prime} ,{\bf q}^{\prime} ,0) \\ \quad =\frac{1}{N(t)}{\rm exp} \left( \frac{{\rm i}}{\hbar }\int _{0}^{t}{\rm d}s\;\left\{ \mathop{\sum }\limits_{\alpha =1}^{2}\left[ {{m}_{\alpha }}{{{\dot{Q}}}_{\alpha }}(s){{{\dot{q}}}_{\alpha }}(s)-{{m}_{\alpha }}\omega _{\alpha }^{2}{{Q}_{\alpha }}(s){{q}_{\alpha }}(s)-{{m}_{\alpha }}Q_{\alpha }^{\prime }{{\gamma }_{\alpha }}(s){{q}_{\alpha }}(s) \right. \right. \right. \\ \qquad -\left. \left. \left. \int _{0}^{s}{\rm d}u\;{{\gamma }_{\alpha }}(s-u){{{\dot{Q}}}_{\alpha }}(u){{q}_{\alpha }}(s) \right]-c(s)\left[ {{Q}_{1}}(s){{q}_{2}}(s)+{{Q}_{2}}(s){{q}_{1}}(s) \right] \right\} \right) \\ \qquad \times {\rm exp} \left\{ -\frac{1}{\hbar }\int _{0}^{t}{\rm d}s\;\int _{0}^{s}{\rm d}u\;\mathop{\sum }\limits_{\alpha =1}^{2}{{K}_{\alpha }}(u-s){{q}_{\alpha }}(s){{q}_{\alpha }}(u) \right\}, \\ \end{array}\end{eqnarray} \tag{ A.37 }$$

where ${\bf Q}=({{Q}_{1}},{{Q}_{2}})$ and ${\bf q}=({{q}_{1}},{{q}_{2}})$. The next step is to determine and solve the equations of motion for the coordinates ${{Q}_{\alpha }}$ and ${{q}_{\alpha }}$. For this, it is necessary to use the Euler–Lagrange equations for the Lagrangian L that appears as the integrand in the imaginary phase of the exponential in equation (A.37). For the coordinates ${{Q}_{\alpha }}$, the equations of motion have the following form

$$\begin{eqnarray}&&\frac{{\rm d}}{{\rm d}s}\left( \frac{\partial L}{\partial {{{\dot{q}}}_{\alpha }}} \right)-\frac{\partial L}{\partial {{q}_{\alpha }}}=0,\end{eqnarray} \tag{ A.38 }$$

whereas for the coordinates ${{q}_{\alpha }}$, the Euler–Lagrange equations are

$$\begin{eqnarray}&&\frac{{\rm d}}{{\rm d}s}\left( \frac{\partial L}{\partial {{{\dot{Q}}}_{\alpha }}} \right)-\frac{\partial L}{\partial {{Q}_{\alpha }}}=0.\end{eqnarray} \tag{ A.39 }$$

By using equation (A.38) in (A.37), the equation of motion for the coordinate ${{Q}_{\alpha }}$ can be written as

$$\begin{eqnarray*}\begin{array}{rcl} {{\ddot{Q}}_{1}}(s)+\omega _{1}^{2}{{Q}_{1}}(s)+\frac{c(s)}{{{m}_{1}}}{{Q}_{2}}(s)+\frac{{\rm d}}{{\rm d}s}\int _{0}^{s}{\rm d}u\;{{\gamma }_{1}}(s-u){{Q}_{1}}(u) & = & 0, \\ \end{array}\end{eqnarray*}$$

$$\begin{eqnarray}\begin{array}{rcl} {{\ddot{Q}}_{2}}(s)+\omega _{2}^{2}{{Q}_{2}}(s)+\frac{c(s)}{{{m}_{2}}}{{Q}_{1}}(s)+\frac{{\rm d}}{{\rm d}s}\int _{0}^{s}{\rm d}u\;{{\gamma }_{2}}(s-u){{Q}_{2}}(u) & = & 0. \\ \end{array}\end{eqnarray} \tag{ A.40 }$$

In order to find the equations of motion for the coordinates ${{q}_{\alpha }}$, consider the following equality

$$\begin{eqnarray}&&\begin{array}{l} Q_{\alpha }^{\prime }\int _{0}^{t}{\rm d}s\;{{\gamma }_{\alpha }}(s){{q}_{\alpha }}(s)+\int _{0}^{t}{\rm d}s\;\int _{0}^{s}{\rm d}u\;{{\gamma }_{\alpha }}(s-u){{{\dot{Q}}}_{\alpha }}(u){{q}_{\alpha }}(s) \\ \quad =q_{\alpha }^{\prime\prime }\int _{0}^{t}{\rm d}s\;{{\gamma }_{\alpha }}(t-s){{Q}_{\alpha }}(s)-\int _{0}^{t}{\rm d}s\;\int _{s}^{t}{\rm d}u\;{{\gamma }_{\alpha }}(u-s){{{\dot{q}}}_{\alpha }}(u){{Q}_{\alpha }}(s). \\ \end{array}\end{eqnarray} \tag{ A.41 }$$

Thus, the propagating function in equation (A.37) can be written as

$$\begin{eqnarray}&&\begin{array}{l} J({\bf Q}^{\prime\prime} ,{\bf q}^{\prime\prime} ,t;{\bf Q}^{\prime} ,{\bf q}^{\prime} ,0) \\ \quad =\frac{1}{N(t)}{\rm exp} \left( \frac{{\rm i}}{\hbar }\int _{0}^{t}{\rm d}s\;\left\{ \mathop{\sum }\limits_{\alpha =1}^{2}\left[ {{m}_{\alpha }}{{{\dot{Q}}}_{\alpha }}(s){{{\dot{q}}}_{\alpha }}(s)-{{m}_{\alpha }}(s)\omega _{\alpha }^{2}{{Q}_{\alpha }}(s){{q}_{\alpha }}(s) \right. \right. \right. \\ \qquad -\left. q_{\alpha }^{\prime\prime }{{\gamma }_{\alpha }}(t-s){{Q}_{\alpha }}(s)+\int _{s}^{t}{\rm d}u\;{{\gamma }_{\alpha }}(u-s){{{\dot{q}}}_{\alpha }}(u){{Q}_{\alpha }}(s) \right] \\ \qquad -\left. \left. c(s)\left[ {{Q}_{1}}(s){{q}_{2}}(s)+{{Q}_{2}}(s){{q}_{1}}(s) \right] \right\} \right) \\ \qquad \times {\rm exp} \left\{ -\frac{1}{\hbar }\int _{0}^{t}{\rm d}s\;\int _{0}^{s}{\rm d}u\;\mathop{\sum }\limits_{\alpha =1}^{2}{{K}_{\alpha }}(u-s){{q}_{\alpha }}(s){{q}_{\alpha }}(u) \right\}. \\ \end{array}\end{eqnarray} \tag{ A.42 }$$

By virtue of the expression (A.39), the equations of motion for the coordinates ${{q}_{\alpha }}$ are given by

$$\begin{eqnarray}\begin{array}{rcl} {{\ddot{q}}_{1}}(s)+\omega _{1}^{2}{{q}_{1}}(s)+\frac{c(s)}{{{m}_{1}}}{{q}_{2}}(s)-\frac{{\rm d}}{{\rm d}s}\int _{s}^{t}{\rm d}u\;{{\gamma }_{1}}(u-s){{q}_{1}}(u) & = & 0, \\ {{\ddot{q}}_{2}}(s)+\omega _{2}^{2}{{q}_{2}}(s)+\frac{c(s)}{{{m}_{2}}}{{q}_{1}}(s)-\frac{{\rm d}}{{\rm d}s}\int _{s}^{t}{\rm d}u\;{{\gamma }_{2}}(u-s){{q}_{2}}(u) & = & 0. \\ \end{array}\end{eqnarray} \tag{ A.43 }$$

The coordinates ${{Q}_{\alpha }}$ and ${{q}_{\alpha }}$ satisfy the boundary conditions

$$\begin{eqnarray}{{Q}_{\alpha }}(s)=\left\{ \begin{array}{ccccccccccccccc} Q_{\alpha }^{\prime }, & s=0, \\ Q_{\alpha }^{\prime\prime }, & s=t, \\ \end{array} \right.\quad {{q}_{\alpha }}(s)=\left\{ \begin{array}{ccccccccccccccc} q_{\alpha }^{\prime }, & s=0, \\ q_{\alpha }^{\prime\prime }, & s=t. \\ \end{array} \right.\end{eqnarray} \tag{ A.44 }$$

Since the equations of motion (A.41) and (A.44) are linear, the solution to them can be written as in equation (13). Once the solutions for the equations of motion are calculated, only one step is left to find the propagating function. The partial integration of the first term in the imaginary phase in the propagating function in (A.37) leads to

$$\begin{eqnarray*}&&\begin{array}{l} J({\bf Q}^{\prime\prime} ,{\bf q}^{\prime\prime} ,t;{\bf Q}^{\prime} ,{\bf q}^{\prime} ,0) \\ \quad =\frac{1}{N(t)}{\rm exp} \left( \frac{{\rm i}}{\hbar }\left\{ q_{1}^{\prime\prime }{{{\dot{Q}}}_{1}}(t)-q_{1}^{\prime }{{{\dot{Q}}}_{1}}(0)+q_{2}^{\prime\prime }{{{\dot{Q}}}_{2}}(t)-q_{2}^{\prime }{{{\dot{Q}}}_{2}}(0) \right. \right. \\ \qquad -{{m}_{1}}\int _{0}^{t}{\rm d}s\;\left[ {{\ddot{Q}}_{1}}(s)+\omega _{1}^{2}{{Q}_{1}}(s)+\frac{c(s)}{{{m}_{1}}}{{Q}_{2}}(s)+Q_{1}^{\prime }{{\gamma }_{1}}(s) \right. \\ \qquad +\left. \int _{0}^{s}{\rm d}u\;{{\gamma }_{1}}(s-u){{{\dot{Q}}}_{1}}(u) \right]{{q}_{1}}(s)-{{m}_{2}}\int _{0}^{t}{\rm d}s\;\left[ {{\ddot{Q}}_{2}}(s)+\omega _{2}^{2}{{Q}_{2}}(s) \right. \\ \qquad +\left. \left. \left. \frac{c(s)}{{{m}_{2}}}{{Q}_{1}}(s)+Q_{2}^{\prime }{{\gamma }_{2}}(s)+\int _{0}^{s}{\rm d}u\;{{\gamma }_{2}}(s-u){{{\dot{Q}}}_{2}}(u) \right]{{q}_{2}}(s) \right\} \right) \\ \qquad \times {\rm exp} \left\{ -\frac{1}{\hbar }\int _{0}^{t}{\rm d}s\;\int _{0}^{s}{\rm d}u\;\mathop{\sum }\limits_{\alpha =1}^{2}{{K}_{\alpha }}(u-s){{q}_{\alpha }}(s){{q}_{\alpha }}(u) \right\}. \\ \end{array}\end{eqnarray*}$$

The last two terms inside the integrals in the imaginary phase of the exponential are the classical equations of motion in (A.41) and (A.43); therefore, their contribution vanishes. Hence, the propagating function takes the compact form

$$\begin{eqnarray}&&J({\bf Q}^{\prime\prime} ,{\bf q}^{\prime\prime} ,t;{\bf Q}^{\prime} ,{\bf q}^{\prime} ,0)=\frac{1}{N(t)}{\rm exp} \left\{ \frac{{\rm i}}{\hbar }\mathop{\sum }\limits_{\alpha =1}^{2}\left[ q_{\alpha }^{\prime\prime }{{{\dot{Q}}}_{\alpha }}(t)-q_{\alpha }^{\prime }{{{\dot{Q}}}_{\alpha }}(0) \right] \right.\end{eqnarray} \tag{ A.45 }$$

$$\begin{eqnarray}&&\left. -\frac{1}{\hbar }\int _{0}^{t}{\rm d}s\;\int _{0}^{s}{\rm d}u\;\mathop{\sum }\limits_{\alpha =1}^{2}{{K}_{\alpha }}(u-s){{q}_{\alpha }}(s){{q}_{\alpha }}(u) \right\}.\end{eqnarray} \tag{ A.46 }$$

By replacing equation (13) into the last expression, the propagating function reads

$$\begin{eqnarray}&&\begin{array}{l} J({\bf Q}^{\prime\prime} ,{\bf q}^{\prime\prime} ,t;{\bf Q}^{\prime} ,{\bf q}^{\prime} ,0) \\ \quad =\frac{1}{N(t)}{\rm exp} \left\{ \frac{{\rm i}}{\hbar }\left[ q_{1}^{\prime\prime }\left( {{{\dot{U}}}_{1}}(t,t)Q_{1}^{\prime }+{{{\dot{U}}}_{2}}(t,t)Q_{1}^{\prime\prime }+{{{\dot{U}}}_{3}}(t,t)Q_{2}^{\prime }+{{{\dot{U}}}_{4}}(t,t)Q_{2}^{\prime\prime } \right) \right. \right. \\ \qquad -q_{1}^{\prime }\left( {{{\dot{U}}}_{1}}(t,0)Q_{1}^{\prime }+{{{\dot{U}}}_{2}}(t,t)Q_{1}^{\prime\prime }+{{{\dot{U}}}_{3}}(t,0)Q_{2}^{\prime }+{{{\dot{U}}}_{4}}(t,0)Q_{2}^{\prime\prime } \right) \\ \qquad +q_{2}^{\prime\prime }\left( {{{\dot{V}}}_{1}}(t,t)Q_{2}^{\prime }+{{{\dot{V}}}_{2}}(t,t)Q_{2}^{\prime\prime }+{{{\dot{V}}}_{3}}(t,t)Q_{1}^{\prime }+{{{\dot{V}}}_{4}}(t,t)Q_{1}^{\prime\prime } \right) \\ \qquad -\left. \left. q_{2}^{\prime }\left( {{{\dot{V}}}_{1}}(t,0)Q_{2}^{\prime }+{{{\dot{V}}}_{2}}(t,0)Q_{2}^{\prime\prime }+{{{\dot{V}}}_{3}}(t,0)Q_{1}^{\prime }+{{{\dot{V}}}_{4}}(t,0)Q_{1}^{\prime\prime } \right) \right] \right\} \\ \qquad \times {\rm exp} \left\{ -\frac{1}{\hbar }\int _{0}^{t}{\rm d}s\;\int _{0}^{s}{\rm d}u\;{{K}_{1}}(u-s)\left[ {{u}_{1}}(t,s)q_{1}^{\prime }+{{u}_{2}}(t,s)q_{1}^{\prime\prime }+{{u}_{3}}(t,s)q_{2}^{\prime } \right. \right. \\ \qquad +\left. \left. {{u}_{4}}(t,s)q_{2}^{\prime\prime } \right]\times \left[ {{u}_{1}}(t,u)q_{1}^{\prime }+{{u}_{2}}(t,u)q_{1}^{\prime\prime }+{{u}_{3}}(t,u)q_{2}^{\prime }+{{u}_{4}}(t,u)q_{2}^{\prime\prime } \right] \right\} \\ \qquad \times {\rm exp} \left\{ -\frac{1}{\hbar }\int _{0}^{t}{\rm d}s\;\int _{0}^{s}{\rm d}u\;{{K}_{2}}(u-s)\left[ {{v}_{1}}(t,s)q_{2}^{\prime }+{{v}_{2}}(t,s)q_{2}^{\prime\prime }+{{v}_{3}}(t,s)q_{1}^{\prime } \right. \right. \\ \qquad +\left. \left. {{v}_{4}}(t,s)q_{1}^{\prime\prime } \right]\times \left[ {{v}_{1}}(t,u)q_{2}^{\prime }+{{v}_{2}}(t,u)q_{2}^{\prime\prime }+{{v}_{3}}(t,u)q_{1}^{\prime }+{{v}_{4}}(t,u)q_{1}^{\prime\prime } \right] \right\}, \\ \end{array}\end{eqnarray} \tag{ A.47 }$$

or alternatively,

$$\begin{eqnarray}&&J({\bf Q}^{\prime\prime} ,{\bf q}^{\prime\prime} ,t;{\bf Q}^{\prime} ,{\bf q}^{\prime} ,0)=\frac{1}{N(t)}{\rm exp} \left\{ \frac{{\rm i}}{\hbar }{{{\boldsymbol{x}} }_{1}}\cdot {\sf A}\cdot {{{\boldsymbol{x}} }_{1}}-\frac{1}{\hbar }{{{\boldsymbol{x}} }_{2}}\cdot {\sf B} \cdot {{{\boldsymbol{x}} }_{2}} \right\},\end{eqnarray} \tag{ A.48 }$$

where ${{{\boldsymbol{x}} }_{1}}=({\bf Q}^{\prime\prime} ,{\bf Q}^{\prime} ,{\bf q}^{\prime\prime} ,{\bf q}^{\prime} )$, ${{{\boldsymbol{x}} }_{2}}=({\bf q}^{\prime\prime} ,{\bf q}^{\prime} )$ and the matrix elements A_ij and B_ij of the matrices ${\sf A}$ and ${\sf B}$ are given by

$$\begin{eqnarray}&&{{A}_{ij}}=\frac{1}{2}\frac{\partial }{\partial x_{1}^{i}}\frac{\partial {{\phi }_{{\rm I}}}}{\partial x_{1}^{j}},\qquad {{B}_{ij}}=\frac{1}{2}\frac{\partial }{\partial x_{2}^{i}}\frac{\partial {{\phi }_{{\rm R}}}}{\partial x_{2}^{j}},\end{eqnarray} \tag{ A.49 }$$

respectively. Here ${{\phi }_{{\rm I}}}$ and ${{\phi }_{{\rm 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

$$\begin{eqnarray}\end{eqnarray} \tag{ A.50 }$$