Heat transport in harmonic oscillator systems with thermal baths: application to optomechanical arrays

Heat conduction in a physical system is a notoriously complex issue to investigate, as the dynamics depend strongly on the interaction between the system constituents as well as on the nature of the environmental baths and their coupling with the system [1–8]. However, low-dimensional systems in contact with different energy or particle baths represent an excellent test-bed for some of the most recent ideas in classical and quantum out-of-equilibrium statistical physics [1]. For example, one can show that a chain of quantum harmonic oscillators in contact with two heat baths at different temperatures exhibits a steady state fluctuation theorem, setting constraints on the entropy production [9], in all respects equivalent to the fluctuation theorem for the corresponding classical case [10]. Furthermore, chains of oscillators have been used as model systems to study heat conduction in solids, in particular to test the validity of Fourier law, according to which the heat current across a material subject to a temperature gradient scales as the inverse of the system size [1]. Motivated by the growing interest in the thermodynamic properties of out-of-equilibrium quantum systems [11–13], we investigate in this paper a prototypical system consisting of a set of quantum harmonic oscillators, each in contact with an independent thermal bath, and all coupled to a common oscillator, which is itself in contact with its own bath. The coupling to the common oscillator effectively mediates an interaction between the different oscillators and baths, which renders the description of the quantum dynamics quite complex in general. While one can envisage many situations that this system may model, our study is specifically motivated by opto- or electromechanical arrays [14–20, 22, 23], in which electromagnetic radiation can affect the motion of mechanically compliant structures, thereby allowing effective transport of phonons between the mechanical elements [25]. In addition to their widespread use for sensing and for communication technologies, opto/electromechanical systems have made great progress towards operation in the quantum regime in the past decade [23]. This has motivated, among other things, their potential application to quantum thermodynamics and the investigation of quantum heat engines, pistons, etc [24–29]. Among their chief virtues is the highly tunable coupling with electromagnetic radiation which can be enhanced with a resonator and which allows for flexible engineering of interactions and readout. Arrays of mechanical oscillators are particularly interesting as long-range interactions between the mechanics can be engineered [20, 21] in a well-controlled fashion and collective phenomena, such as self-oscillations, phonon lasing or synchronization, can occur [15, 62, 31–36].

In this work we investigate phonon transport in an ensemble of identical oscillators in contact with independent thermal baths with (possibly) different temperatures and coupled to a common oscillator—a single electromagnetic field mode in an opto/electromechanical array setting—as illustrated schematically in figure 1. We shall follow the standard treatment to adiabatically-eliminate this common oscillator. Under these conditions, this gives rise to a single shared bath for the rest of the system, which can be mathematically recast without further approximation as a number of individual correlated baths, one connected to each oscillator. Our calculation is valid in the regime typical of optomechanical experiments, where the approximations made (e.g., adiabatic elimination of the optical mode) work very well, but let us remark that the problem could alternatively be tackled by means of the scattering matrix theory or the Landauer approach [1]. A full comparison between our results and the exact results obtained using these methods is worth pursuing but lies outside the scope of the present article.

Zoom In Zoom Out Reset image size

We start by deriving a general expression, equation (9), for the heat flow through the individual elements in steady state when the couplings to the common mode are weak and this mode can be adiabatically eliminated. After discussing and solving the special case where all the baths are at the same temperature (section 3), we consider the general problem consisting of a two-mechanical-resonator array coupled to one electromagnetic field and thermal baths with arbitrary temperatures. We show that, after the adiabatic elimination of the field, it is equivalent to a generic two-oscillator system with an effective mutual linear coupling, an effective common bath and two independent baths (sections 4.1 and 4.2). We solve this generic problem for thermal Markovian baths and derive expressions, equations (39) and(40), for the steady state occupation and heat flows of the mechanics. In section 4.3, we discuss the results in various parameter regimes which could be realized through a suitable engineering of the optomechanical interaction and give an indication of the various systems that could be used for investigating the effects we explore. Our results show the possibility for engineering heat flow and the Fourier law in arrays of harmonic oscillators possessing only indirect coupling, and we illustrate this by proposing a practical application in optomechanical arrays. Systems with several oscillators, where the interaction between the different constituents can be easily tuned, and the local temperature can be precisely controlled are highly desirable given the current interest in out-of equilibrium physics, and so the set-up we propose can represent an advancement with respect to the experimental set-ups which have been recently used to study the entropy production in systems in contact with only two heat baths [37, 38]. Likewise, our system can represent a test-bed to extend some of the concepts of stochastic thermodynamics [39] to the quantum case, allowing one to measure, e.g., the heat current, and thus the entropy production in a quantum out-of-equilibrium system. Finally, we conclude in section 5 by surveying our results and putting them in the context of possible future work.

We consider a system composed of N identical harmonic oscillators, mutually uncoupled, but all linearly coupled to a common harmonic oscillator, and in contact with independent thermal baths. We denote by ϱ the density matrix for the $(N+1)$-partite system. The master equation for ϱ can be written ($\hbar =1$) :

$$\begin{eqnarray}&&\dot{\varrho }=-{\rm i}[\hat{H},\varrho ]+\mathop{\sum }\limits_{j=1}^{N}{{\mathcal{L}}_{j}}\varrho +{{\mathcal{L}}_{{\rm a}}}\varrho ,\end{eqnarray} \tag{ 1 }$$

where $\hat{H}$ is the Hamiltonian governing the evolution,

$$\begin{eqnarray}&&{{\mathcal{L}}_{j}}\varrho ={{\gamma }_{j}}({{n}_{j}}+1)D[{{\hat{b}}_{j}}]\varrho +{{\gamma }_{j}}{{n}_{j}}D[\hat{b}_{j}^{\dagger }]\varrho ,\end{eqnarray} \tag{ 2 }$$

with

$$\begin{eqnarray}&&D[\hat{b}]\varrho =2\hat{b}\varrho {{\hat{b}}^{\dagger }}-{{\hat{b}}^{\dagger }}\hat{b}\varrho -\varrho {{\hat{b}}^{\dagger }}\hat{b}\end{eqnarray} \tag{ 3 }$$

and ${{\hat{b}}_{j}}$ the annihilation operator of the jth oscillator ($1\leqslant j\leqslant N$), ${{\gamma }_{j}}$ the coupling rate of the jth oscillator to its bath, whose mean occupation number is n_j, and

$$\begin{eqnarray}&&{{\mathcal{L}}_{{\rm a}}}\varrho =\kappa D[\hat{a}]\varrho ,\end{eqnarray} \tag{ 4 }$$

with $\hat{a}$ the annihilation operator of the final harmonic oscillator, and κ its coupling constant to its own bath, assumed at zero-temperature. We assume that $\hat{H}$ has the following form [25]:

$$\begin{eqnarray}&&\hat{H}=\Omega {{\hat{a}}^{\dagger }}\hat{a}+\mathop{\sum }\limits_{j=1}^{N}\omega \hat{b}_{j}^{\dagger }{{\hat{b}}_{j}}+\mathop{\sum }\limits_{j=1}^{N}{{g}_{j}}\left( \hat{a}+{{{\hat{a}}}^{\dagger }} \right)\left( {{{\hat{b}}}_{j}}+\hat{b}_{j}^{\dagger } \right),\end{eqnarray} \tag{ 5 }$$

where ω and Ω are the oscillators' frequencies and g_j, assumed real, represents the coupling strength of the jth oscillator to the final one. If g_j is small compared to the other frequencies of the problem, $\hat{a}$ can be adiabatically eliminated from the dynamics⁴ . For the interested reader, we reproduce this elimination procedure in appendix A.

Let us indicate the reduced, N-partite, density matrix that results from this elimination process with ρ. The heat flow into or out of the lth element⁵ is given by [5]

$$\begin{eqnarray}&&{{J}_{l}}={\rm Tr}\left( \hat{H}{{\mathcal{L}}_{l}}\varrho \right).\end{eqnarray} \tag{ 6 }$$

We now make use of our adiabatic elimination procedure to write, in steady state and to lowest order in the coupling constants g_j, $\varrho ={{\rho }_{{\rm ss}}}\;\otimes \;{{\rho }_{{\rm a}}}$, where ${{\rho }_{{\rm a}}}$ is the steady-state density matrix for the $(N+1){\rm st}$ harmonic oscillator. Let us first take the trace with respect to this mode⁶ :

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{J}_{l}} & = & {\rm T}{{{\rm r}}_{\neg {\rm a}}}\left( {\rm T}{{{\rm r}}_{{\rm a}}}\left\{ \left[ \Omega {{{\hat{a}}}^{\dagger }}\hat{a}+\mathop{\sum }\limits_{j=1}^{N}\omega \hat{b}_{j}^{\dagger }{{{\hat{b}}}_{j}}+\mathop{\sum }\limits_{j=1}^{N}{{g}_{j}}\left( \hat{a}+{{{\hat{a}}}^{\dagger }} \right)\left( {{{\hat{b}}}_{j}}+\hat{b}_{j}^{\dagger } \right) \right] \right. \right. \\ {} & {} & \times \left. \left. {{\mathcal{L}}_{l}}{{\rho }_{{\rm ss}}}\;\otimes \;{{\rho }_{{\rm a}}} \right\} \right) \\ {} & = & {\rm T}{{{\rm r}}_{\neg {\rm a}}}\left\{ \left[ \Omega \langle {{{\hat{a}}}^{\dagger }}\hat{a}\rangle +\mathop{\sum }\limits_{j=1}^{N}\omega \hat{b}_{j}^{\dagger }{{{\hat{b}}}_{j}}+\mathop{\sum }\limits_{j=1}^{N}{{g}_{j}}\langle \hat{a}+{{{\hat{a}}}^{\dagger }}\rangle \left( {{{\hat{b}}}_{j}}+\hat{b}_{j}^{\dagger } \right) \right]{{\mathcal{L}}_{l}}{{\rho }_{{\rm ss}}} \right\}. \\ \end{array}\end{eqnarray} \tag{ 7 }$$

Following the usual methodology, and under the assumption that the oscillator described by $\hat{a}$ is driven into a coherent state, the state described by ${{\rho }_{{\rm a}}}$ has been shifted to describe a zero-mean vacuum state, so that $\langle \hat{a}+{{\hat{a}}^{\dagger }}\rangle =0$ and $\langle {{\hat{a}}^{\dagger }}\hat{a}\rangle =0$. The adiabatic elimination used to derive the above expression assumes a weak coupling between the two subsystems. Consistently with this approximation, we keep our results at lowest order in the coupling strength and neglect the effects of correlations between mode 'a' and the other modes. We then obtain

$$\begin{eqnarray}&&{{J}_{l}}={\rm Tr}\left\{ \omega \mathop{\sum }\limits_{j=1}^{N}\hat{b}_{j}^{\dagger }{{{\hat{b}}}_{j}}{{\mathcal{L}}_{l}}{{\rho }_{{\rm ss}}} \right\},\end{eqnarray} \tag{ 8 }$$

where we drop the subscript from the trace because there is no longer any ambiguity. By exploiting the bosonic commutation relations and the cyclic property of the trace, we obtain an explicit expression for the steady-state heat flow:

$$\begin{eqnarray}&&{{J}_{l}}=2\omega {{\gamma }_{l}}[{{n}_{l}}-{\rm Tr}(\hat{b}_{l}^{\dagger }{{\hat{b}}_{l}}{{\rho }_{{\rm ss}}})]=2\omega {{\gamma }_{l}}\left( {{n}_{l}}-\langle \hat{b}_{l}^{\dagger }{{{\hat{b}}}_{l}}\rangle \right).\end{eqnarray} \tag{ 9 }$$

The work in this section extends the treatments of [5] and [40] to the case where the N oscillators do not interact directly, but via coupling to a common, adiabatically eliminated, mode. To proceed further, one needs an expression for the average occupation of the mechanical elements. As we shall now illustrate, there are various situations under which we can calculate this quantity explicitly.

It is possible to obtain an analytic expression for the heat flow in the case in which the different thermal baths to which the elements are connected are identical, i.e., they are characterized by one single coupling constant ${{\gamma }_{j}}=\gamma$ and occupation number ${{n}_{j}}=n$. To continue, it is convenient to introduce a normal-mode basis. Define $g\sqrt{\sum _{j=1}^{N}g_{j}^{2}}\ne 0$, ${{{\boldsymbol{g}} }^{(1)}}={{({{g}_{1}}/g\ {{g}_{2}}/g\ ...\ {{g}_{N}}/g)}^{{\rm T}}}$, and ${{{\boldsymbol{g}} }^{(l)}}$ ($2\leqslant l\leqslant N$) such that the set $\{{{{\boldsymbol{g}} }^{(l)}}\}$ forms an orthonormal basis⁷ . We use this 'collective' basis to define a new set of N normalized harmonic oscillator annihilation operators ${{\tilde{b}}_{j}}$ ($1\leqslant j\leqslant N$):

$$\begin{eqnarray}&&\left( \begin{array}{ccccccccccccccc} {{\tilde{b}}_{1}} \\ {{\tilde{b}}_{2}} \\ \vdots \\ {{\tilde{b}}_{N}} \\ \end{array} \right)=\left[ \begin{array}{ccccccccccccccc} {{{\boldsymbol{g}} }^{(1)}}^{{\rm T}} \\ {{{\boldsymbol{g}} }^{(2)}}^{{\rm T}} \\ \vdots \\ {{{\boldsymbol{g}} }^{(N)}}^{{\rm T}} \\ \end{array} \right]\cdot \left( \begin{array}{ccccccccccccccc} {{{\hat{b}}}_{1}} \\ {{{\hat{b}}}_{2}} \\ \vdots \\ {{{\hat{b}}}_{N}} \\ \end{array} \right)={\boldsymbol{G}} \cdot \left( \begin{array}{ccccccccccccccc} {{{\hat{b}}}_{1}} \\ {{{\hat{b}}}_{2}} \\ \vdots \\ {{{\hat{b}}}_{N}} \\ \end{array} \right),\end{eqnarray} \tag{ 10 }$$

whereby

$$\begin{eqnarray}&&\dot{\varrho }=-{\rm i}[\hat{H},\varrho ]+\mathop{\sum }\limits_{j=1}^{N}{{\tilde{\mathcal{L}}}_{j}}\varrho +{{\mathcal{L}}_{{\rm a}}}\varrho ,\end{eqnarray} \tag{ 11 }$$

where each ${{\tilde{\mathcal{L}}}_{j}}$ is defined analogously to ${{\mathcal{L}}_{j}}$ but with ${{\hat{b}}_{j}}$ replaced by ${{\tilde{b}}_{j}}$. The Hamiltonian can be expressed as the sum of Hamiltonians for $N-1$ uncoupled and two linearly-coupled oscillators:

$$\begin{eqnarray}&&\hat{H}=\Omega {{\hat{a}}^{\dagger }}\hat{a}+\mathop{\sum }\limits_{j=1}^{N}\omega \tilde{b}_{j}^{\dagger }{{\tilde{b}}_{j}}+g\left( \hat{a}+{{{\hat{a}}}^{\dagger }} \right)\left( {{\tilde{b}}_{1}}+\tilde{b}_{1}^{\dagger } \right).\end{eqnarray} \tag{ 12 }$$

Note that it is only because the thermal baths all have the same temperature that the Liouvillians ${{\tilde{L}}_{j}}$ are diagonal in the new basis. As shown in appendix A, the adiabatic elimination of $\hat{a}$ and application of the rotating-wave approximation leads to a shift of the oscillator ${{\tilde{b}}^{1}}$'s frequency and a modification of its coupling to the baths. In the single-oscillator case the adiabatic elimination leads to the appearance of an effective thermal bath. In the new basis, this means that we obtain an effective master equation for the reduced oscillator-only subsystem

$$\begin{eqnarray}&&\dot{\rho }=-{\rm i}[{{\hat{H}}_{{\rm eff}}},\rho ]+{{\tilde{\mathcal{L}}}^{\prime }}\rho +\mathop{\sum }\limits_{j=2}^{N}{{\tilde{\mathcal{L}}}_{j}}\rho ,\end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray}&&{{\hat{H}}_{{\rm eff}}}={{\omega }^{\prime }}\tilde{b}_{1}^{\dagger }{{\tilde{b}}_{1}}+\mathop{\sum }\limits_{j=2}^{N}\omega \tilde{b}_{j}^{\dagger }{{\tilde{b}}_{j}},\end{eqnarray} \tag{ 14 }$$

and

$$\begin{eqnarray}&&{{\tilde{\mathcal{L}}}^{\prime }}\rho ={{\gamma }^{\prime }}\left[ ({{n}^{\prime }}+1)\left( 2{{\tilde{b}}_{1}}\rho \tilde{b}_{1}^{\dagger }-\tilde{b}_{1}^{\dagger }{{\tilde{b}}_{1}}\rho -\rho \tilde{b}_{1}^{\dagger }{{\tilde{b}}_{1}} \right)+{{n}^{\prime }}\left( 2\tilde{b}_{1}^{\dagger }\rho {{\tilde{b}}_{1}}-{{\tilde{b}}_{1}}\tilde{b}_{1}^{\dagger }\rho -\rho {{\tilde{b}}_{1}}\tilde{b}_{1}^{\dagger } \right) \right].\end{eqnarray} \tag{ 15 }$$

The reduced master equation (13) has a steady-state solution ${{\rho }_{{\rm ss}}}$ given by a tensor product of N thermal states ${{\rho }_{j}}$. We define ${{n}^{\prime }}=\langle \tilde{b}_{1}^{\dagger }{{\tilde{b}}_{1}}\rangle$ to be the occupation number of the state for j = 1, and $n=\langle \tilde{b}_{j}^{\dagger }{{\tilde{b}}_{j}}\rangle$ for $j\geqslant 2$. Writing the heat flow in the normal-mode basis yields

$$\begin{eqnarray}&&{{J}_{l}}=2\omega \gamma \left( n-\mathop{\sum }\limits_{j,k=1}^{N}{{{\boldsymbol{G}} }_{jl}}{{{\boldsymbol{G}} }_{kl}}\langle \tilde{b}_{j}^{\dagger }{{\tilde{b}}_{k}}\rangle \right).\end{eqnarray} \tag{ 16 }$$

Because of the uncorrelated nature of the steady state, we can immediately write that $\langle \tilde{b}_{j}^{\dagger }{{\tilde{b}}_{k}}\rangle =0$ if $j\ne k$, i.e.

$$\begin{eqnarray}&&{{J}_{l}}=2\omega \gamma \mathop{\sum }\limits_{j=1}^{N}{\boldsymbol{G}} _{jl}^{2}\left( n-\langle \tilde{b}_{j}^{\dagger }{{\tilde{b}}_{j}}\rangle \right),\end{eqnarray} \tag{ 17 }$$

where we have also taken the sum out of the parentheses by exploiting the properties of the orthonormal matrix ${\boldsymbol{G}}$. By the very nature of the steady state, however,

$$\begin{eqnarray}&&\langle \tilde{b}_{j}^{\dagger }{{\tilde{b}}_{j}}\rangle =n+({{n}^{\prime }}-n){{\delta }_{j,1}},\end{eqnarray} \tag{ 18 }$$

so that

$$\begin{eqnarray}&&{{J}_{l}}=2\omega \gamma {\boldsymbol{G}} _{1\,l}^{2}(n-{{n}^{\prime }})=2\omega \gamma (n-{{n}^{\prime }})g_{l}^{2}/{{g}^{2}}.\end{eqnarray} \tag{ 19 }$$

Therefore, the heat flowing into or out of the mechanical subsystem is simply

$$\begin{eqnarray}&&{{J}_{{\rm m}}}\mathop{\sum }\limits_{l=1}^{N}{{J}_{l}}=2\omega \gamma (n-{{n}^{\prime }}),\end{eqnarray} \tag{ 20 }$$

which is nonzero because we have traced out the $(N+1){\rm st}$ oscillator. Let us note that $n-{{n}^{\prime }}\to 0$ when $g\to 0$. The heat flowing into or out of this oscillator must therefore be

$$\begin{eqnarray}&&{{J}_{{\rm c}}}-\mathop{\sum }\limits_{l=1}^{N}{{J}_{l}}=2\omega \gamma ({{n}^{\prime }}-n),\end{eqnarray} \tag{ 21 }$$

to maintain balance, i.e., ${{J}_{{\rm m}}}+{{J}_{{\rm c}}}=0$.

4.1. Two-element optomechanical array

As an application, we consider a system composed of two mechanical oscillators in which each identical, independent oscillator is dispersively coupled by radiation pressure to the same cavity field mode. The mechanical oscillators have identical frequency ω and equal damping rate γ into two independent Markovian thermal baths held at possibly different temperatures, yielding mean thermal occupation numbers n₁ and n₂ for the mechanical modes in absence of the field. We assume operation in the linearized regime for the optomechanical interaction [23] in which the number of intracavity photons is large. Without loss of generality we also consider a situation in which the cavity field couples to the centre-of-mass motion of the pair of mechanical oscillators, leaving the relative mode of motion uncoupled⁸ . Under these conditions the Hamiltonian and Liouvillian of the system are given by

$$\begin{eqnarray}&&H=-\Delta {{\hat{a}}^{\dagger }}\hat{a}+\omega (\hat{b}_{1}^{\dagger }{{\hat{b}}_{1}}+\hat{b}_{2}^{\dagger }{{\hat{b}}_{2}})+g(\hat{a}+{{\hat{a}}^{\dagger }})\left[ ({{{\hat{b}}}_{1}}+\hat{b}_{1}^{\dagger })+({{{\hat{b}}}_{2}}+\hat{b}_{2}^{\dagger }) \right]/\sqrt{2},\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&\mathcal{L}={{\mathcal{L}}_{1}}+{{\mathcal{L}}_{2}}+{{\mathcal{L}}_{a}}.\end{eqnarray} \tag{ 23 }$$

We have defined the cavity field detuning $\Delta ={{\Omega }_{{\rm L}}}-\Omega$, where ${{\Omega }_{{\rm L}}}$ is the frequency of the driving field, assumed to be monochromatic, and have transformed our system to a frame rotating at the frequency ${{\Omega }_{{\rm L}}}$. The enhanced optomechanical coupling rate g and the Liouvillians ${{\mathcal{L}}_{1(2)}}$ and ${{\mathcal{L}}_{{\rm a}}}$ are defined as in equations (2) and (4), respectively. Introducing the relative and centre-of-mass modes

$$\begin{eqnarray}&&{{\hat{b}}_{{\rm r}}}=({{\hat{b}}_{1}}-{{\hat{b}}_{2}})/\sqrt{2},\qquad {{\hat{b}}_{{\rm c}}}=({{\hat{b}}_{1}}+{{\hat{b}}_{2}})/\sqrt{2},\end{eqnarray} \tag{ 24 }$$

these can be recast as

$$\begin{eqnarray}&&H=-\Delta {{\hat{a}}^{\dagger }}\hat{a}+\omega (\hat{b}_{{\rm r}}^{\dagger }{{\hat{b}}_{{\rm r}}}+\hat{b}_{{\rm c}}^{\dagger }{{\hat{b}}_{{\rm c}}})+g(\hat{a}+{{\hat{a}}^{\dagger }})({{\hat{b}}_{{\rm c}}}+\hat{b}_{{\rm c}}^{\dagger }),\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&\mathcal{L}={{\mathcal{L}}_{{\rm r}}}+{{\mathcal{L}}_{{\rm c}}}+{{\mathcal{L}}_{{\rm r},{\rm c}}}+{{\mathcal{L}}_{{\rm a}}},\end{eqnarray} \tag{ 26 }$$

with ${{\mathcal{L}}_{{\rm r},{\rm c}}}$ a Liouvillian including the correlations between the rotated baths.

Adiabatically eliminating the cavity field in the weak optomechanical coupling regime leads to an effective mechanical frequency ${{\omega }^{\prime }}=\omega +\Lambda$, due to the optical spring effect, as well as an effective damping rate ${{\gamma }^{\prime }}=\gamma +\bar{\gamma }$, for the centre-of-mass mode (see appendix A), where

$$\begin{eqnarray}&&\Lambda =\frac{2{{g}^{2}}\Delta ({{\Delta }^{2}}-{{\omega }^{2}}+{{\kappa }^{2}})}{[{{(\Delta +\omega )}^{2}}+{{\kappa }^{2}}][{{(\Delta -\omega )}^{2}}+{{\kappa }^{2}}]},\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&\bar{\gamma }={{\beta }_{+}}-{{\beta }_{-}}=\frac{{{G}^{2}}\kappa }{{{(\Delta +\omega )}^{2}}+{{\kappa }^{2}}}-\frac{{{G}^{2}}\kappa }{{{(\Delta -\omega )}^{2}}+{{\kappa }^{2}}}.\end{eqnarray} \tag{ 28 }$$

The field fluctuations also give rise to an additional effective bath, coupling to the centre-of-mass mode only, whose Liouvillian has the form

$$\begin{eqnarray}&&\bar{\mathcal{L}}={{\beta }_{+}}D[{{\hat{b}}_{{\rm c}}}]+{{\beta }_{-}}D[\hat{b}_{{\rm c}}^{\dagger }].\end{eqnarray} \tag{ 29 }$$

Assuming a red-detuned cavity field ($\Delta \lt 0$) for which ${{\beta }_{+}}\gt {{\beta }_{-}}$, this term describes the coupling with a thermal bath with coupling rate and occupation number defined by

$$\begin{eqnarray}&&{{\beta }_{+}}=\bar{\gamma }(\bar{n}+1),\quad {{\beta }_{-}}=\bar{\gamma }\bar{n}\quad \Rightarrow \quad \bar{\gamma }={{\beta }_{+}}-{{\beta }_{-}},\quad \bar{n}=\frac{{{\beta }_{-}}}{{{\beta }_{+}}-{{\beta }_{-}}}.\end{eqnarray} \tag{ 30 }$$

The effective Hamiltonian and Liouvillian after adiabatic elimination then read

$$\begin{eqnarray}&&{{H}_{{\rm eff}}}=(\omega +\Lambda )\hat{b}_{{\rm c}}^{\dagger }{{\hat{b}}_{{\rm c}}}+\omega \hat{b}_{{\rm r}}^{\dagger }{{\hat{b}}_{{\rm r}}},\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&{{\mathcal{L}}_{{\rm eff}}}={{\mathcal{L}}_{{\rm r}}}+{{\mathcal{L}}_{{\rm c}}}+{{\mathcal{L}}_{{\rm r},{\rm c}}}+\bar{\mathcal{L}},\end{eqnarray} \tag{ 32 }$$

or, going back to the bare mechanical basis,

$$\begin{eqnarray}&&{{H}_{{\rm eff}}}=(\omega +\Lambda /2)(\hat{b}_{1}^{\dagger }{{\hat{b}}_{1}}+\hat{b}_{2}^{\dagger }{{\hat{b}}_{2}})+(\Lambda /2)(\hat{b}_{1}^{\dagger }{{\hat{b}}_{2}}+{{\hat{b}}_{1}}\hat{b}_{2}^{\dagger }),\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{{\mathcal{L}}_{{\rm eff}}}={{\mathcal{L}}_{1}}+{{\mathcal{L}}_{2}}+\bar{\mathcal{L}}.\end{eqnarray} \tag{ 34 }$$

In the regime of interest, $\Lambda \ll \omega$, so that we may approximate

$$\begin{eqnarray}&&{{H}_{{\rm eff}}}=\omega (\hat{b}_{1}^{\dagger }{{\hat{b}}_{1}}+\hat{b}_{2}^{\dagger }{{\hat{b}}_{2}})+(\Lambda /2)(\hat{b}_{1}^{\dagger }{{\hat{b}}_{2}}+{{\hat{b}}_{1}}\hat{b}_{2}^{\dagger }).\end{eqnarray} \tag{ 35 }$$

These indeed describe the dynamics of a pair of oscillators, mutually coupled with a strength Λ and in contact with two independent baths and a common bath described by the Liouvillians ${{\mathcal{L}}_{1,2}}$ and $\bar{\mathcal{L}}$, respectively. The master equation that emerges from this model is a special case of a more general master equation; in the next section we shall describe this more general situation and derive the steady state of the system.

4.2. Two oscillators with independent and common thermal baths

A general treatment of the dynamics of a pair of coupled oscillators in contact with either independent baths or a common bath can be found in [4]. Motivated by the results of the previous section, we consider the same Hamiltonian as in [4]:

$$\begin{eqnarray}&&{{\hat{H}}_{{\rm eff}}}=\frac{\omega }{2}\mathop{\sum }\limits_{j=1}^{N}\left( \hat{p}_{j}^{2}+\hat{x}_{j}^{2} \right)+\Lambda {{\hat{x}}_{1}}{{\hat{x}}_{2}},\end{eqnarray} \tag{ 36 }$$

but focus on the case of two identical oscillators (equal mass and frequency) and each in simultaneous contact with both an independent and a common bath; we have defined ${{\hat{x}}_{j}}=({{\hat{b}}_{j}}+\hat{b}_{j}^{\dagger })/\sqrt{2}$ and ${{\hat{p}}_{j}}=({{\hat{b}}_{j}}-\hat{b}_{j}^{\dagger })/({\rm i}\sqrt{2})$ (j = 1, 2). We assume the baths to be Markovian, the independent bath temperatures being given by n_j, and the common bath temperature by $\bar{n}$. In order to obtain an effective master equation in the limit $\gamma ,\bar{\gamma }\ll \Lambda \ll \omega$, we employ these assumptions to follow the procedure of [4], which we will not detail explicitly but finally yields

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {\dot{\rho }} & = & -{\rm i}[{{{\hat{H}}}_{{\rm eff}}},\rho ]-\gamma \mathop{\sum }\limits_{j=1}^{2}\left( {\rm i}[{{{\hat{x}}}_{j}},\{{{{\hat{p}}}_{j}},\rho \}]+(2{{n}_{j}}+1)[{{{\hat{x}}}_{j}},[{{{\hat{x}}}_{j}},\rho ]] \right) \\ {} & {} & -\frac{{\bar{\gamma }}}{2}\mathop{\sum }\limits_{j=1}^{2}\left( {\rm i}[{{{\hat{x}}}_{+}},\{{{{\hat{p}}}_{j}},\rho \}]+(2\bar{n}+1)[{{{\hat{x}}}_{+}},[{{{\hat{x}}}_{j}},\rho ]] \right), \\ \end{array}\end{eqnarray} \tag{ 37 }$$

where ${{\hat{x}}_{+}}={{\hat{x}}_{1}}+{{\hat{x}}_{2}}$. Under the rotating-wave approximation, where terms of the form ${{\hat{b}}_{1}}{{\hat{b}}_{2}}$ and $\hat{b}_{1}^{\dagger }\hat{b}_{2}^{\dagger }$ are neglected based on the approximation that their effects average out over the longer time-scales of relevance to the problem, the model described by this master equation reduces to the optomechanical model derived in section 4.1. In the remainder of this section we show that, under the conditions outlined below, we can solve this more general model explicitly and apply this solution to the situation in section 4.1.

From this master equation, and exploiting the standard commutation relations, we derive a closed system of equations describing the temporal evolution of the sixteen second-order moments, which is shown in appendix B. Under the assumption that the initial state is Gaussian, the dynamics described by these equations preserves the Gaussian nature of the state at all time. From these equations it is straightforward to compute the steady state occupation numbers of both oscillators,

$$\begin{eqnarray}&&n_{j}^{\prime }=\langle \hat{b}_{j}^{\dagger }{{\hat{b}}_{j}}\rangle =\left[ \langle \hat{x}_{j}^{2}\rangle +\langle \hat{p}_{j}^{2}\rangle -1 \right]/2,\end{eqnarray} \tag{ 38 }$$

and thereby calculate the heat flow in this system from (9). One gets

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} n_{j}^{\prime } & = & \frac{2\gamma {{n}_{j}}+\bar{\gamma }\bar{n}}{2\gamma +\bar{\gamma }}+\frac{{{{\bar{\gamma }}}^{2}}}{2(2\gamma +\bar{\gamma })(\gamma +\bar{\gamma })}\left( \frac{{{n}_{1}}+{{n}_{2}}}{2}-\bar{n} \right) \\ {} & {} & +{{(-1)}^{j}}\frac{2\gamma {{\Lambda }^{2}}}{(2\gamma +\bar{\gamma })[{{(2\gamma +\bar{\gamma })}^{2}}+{{\Lambda }^{2}}]}\frac{{{n}_{1}}-{{n}_{2}}}{2}\qquad (j=1,2), \\ \end{array}\end{eqnarray} \tag{ 39 }$$

which leads to the following expression for the steady-state heat flows

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{J}_{j}} & = & 4\omega \gamma \left\{ \frac{{\bar{\gamma }}}{2\gamma +\bar{\gamma }}({{n}_{j}}-\bar{n})+\frac{{{{\bar{\gamma }}}^{2}}}{2(2\gamma +\bar{\gamma })(\gamma +\bar{\gamma })}\left( \bar{n}-\frac{{{n}_{1}}+{{n}_{2}}}{2} \right) \right. \\ {} & {} & -\left. {{(-1)}^{j}}\frac{\gamma {{\Lambda }^{2}}}{(2\gamma +\bar{\gamma })[{{(2\gamma +\bar{\gamma })}^{2}}+{{\Lambda }^{2}}]}({{n}_{1}}-{{n}_{2}}) \right\}\quad (j=1,2). \\ \end{array}\end{eqnarray} \tag{ 40 }$$

These expressions bring together the work of [4] and our generic expression for the heat flow, equation (9). It expresses the heat flow through an array of two oscillators that are coupled not only to independent baths but also to a third, 'correlated' or common, bath, which sets up an effective interaction between the two oscillators, altering the nature of the heat transport. In the next section, we discuss this interplay between the independent and common baths in more detail.

4.3. Discussion

4.3.1. Two-oscillator case

It is interesting to look at these expressions in various regimes of interest. The regime of large mutual coupling between the oscillators corresponds to the large optical spring regime, which can be achieved in the bad-cavity limit of optomechanics, $\kappa \gg \omega$. In this regime, and for $\Delta \sim -\kappa$, equations (27) and (28) give $|\Lambda |\sim {{g}^{2}}/\kappa \gg \bar{\gamma }\sim 2{{g}^{2}}\omega /{{\kappa }^{2}}$. In contrast, the regime of large coupling to the effective common bath can be achieved in the resolved sideband regime of optomechanics, $\omega \gg \kappa$. For $\Delta =-\omega$, one obtains $|\Lambda |\sim {{g}^{2}}/(2\omega )\ll \bar{\gamma }\sim {{g}^{2}}/\kappa$.

• (i)  
  In the large coupling regime $\Lambda \gg \gamma ,\bar{\gamma }$, one has

  $$\begin{eqnarray}&&n_{j}^{\prime }=\frac{2\gamma +\bar{\gamma }}{2(\gamma +\bar{\gamma })}\frac{{{n}_{1}}+{{n}_{2}}}{2}+\frac{{\bar{\gamma }}}{2(\gamma +\bar{\gamma })}\bar{n}\qquad (j=1,2)\end{eqnarray} \tag{ 41 }$$

  which we can understand in two different limits. (a) When the coupling to the independent baths is larger than that to the common bath ($\gamma \gg \bar{\gamma }$), the mean independent bath temperature $({{n}_{1}}+{{n}_{2}})/2$. (b) When the damping into the common bath dominates over the mutual coupling, the whole system equilibrates at the mean of the common and independent bath temperatures

  $$\begin{eqnarray}&&n_{j}^{\prime }=\frac{{\bar{n}}}{2}+\frac{{{n}_{1}}+{{n}_{2}}}{4}\qquad (j=1,2).\end{eqnarray} \tag{ 42 }$$
• (ii)  
  If the independent baths are held at the same temperature n, the term due to the mutual coupling vanishes and one gets

  $$\begin{eqnarray}&&n_{j}^{\prime }=\frac{2\gamma +\bar{\gamma }}{2(\gamma +\bar{\gamma })}n+\frac{{\bar{\gamma }}}{2(\gamma +\bar{\gamma })}\bar{n},\quad {{J}_{j}}=\omega \gamma \frac{{\bar{\gamma }}}{\gamma +\bar{\gamma }}(n-\bar{n})\quad (j=1,2).\end{eqnarray} \tag{ 43 }$$

  This is consistent with the limit of radiation pressure cooling in optomechanics [41–44], which predicts that the (centre-of-mass) mode coupled to the field is cooled down to $\bar{n}$, $\bar{\gamma }\gg \gamma$. Since the uncoupled (relative motion) mode's occupancy remains n, this means that, in the bare basis, $n_{j}^{\prime }\to (n+\bar{n})/2$ when $\bar{\gamma }\gg \gamma$.

As an illustration, figure 2 shows the individual heat flows J₁ and J₂, as well as the total heat flow through the mechanics ${{J}_{1}}+{{J}_{2}}$, as a function of $\Lambda /\gamma$ and $\bar{\gamma }/\gamma$, in the case where the first oscillator is coupled to a higher temperature bath than the second, and for a common bath at zero temperature. The heat flow from the hotter oscillator increases both with Λ and $\bar{\gamma }$, as both the common bath and the other oscillator's bath have a lower temperature. The heat flow from the colder oscillator, however, becomes negative for a large mutual coupling, as this coupling tends to equalize the temperature of both oscillators. Moreover, depending on the temperature difference between the independent baths, the heat flow from the cold oscillator can be seen to either increase (figure 2(b)) or decrease (figure 2(e)) with the coupling with the common bath, as the cold oscillator gets either cooled or heated by the combined action of the cold common bath and the other hot oscillator. The total heat flow through the mechanics is in contrast independent of the mutual coupling and steadily increases with the coupling to the common bath, as the optical field globally takes away heat from the mechanics.

Zoom In Zoom Out Reset image size

4.3.2. N-oscillator case

In the case of N baths all at equal temperature, one can use this formalism together with the ideas developed in section 3 to find

$$\begin{eqnarray}&&n_{j}^{\prime }=n+\frac{g_{j}^{2}}{{{g}^{2}}}\frac{{\bar{\gamma }}}{\gamma +\bar{\gamma }}(\bar{n}-n),\quad {{J}_{j}}=2\omega \gamma \frac{g_{j}^{2}}{{{g}^{2}}}\frac{{\bar{\gamma }}}{\gamma +\bar{\gamma }}(n-\bar{n})\quad (j=1-N).\end{eqnarray} \tag{ 44 }$$

This means that the heat flow through the jth mechanical element is proportional to the temperature difference between the independent thermal baths and the field bath, weighted by the branching ratio of the damping rates $\eta =\bar{\gamma }/(\gamma +\bar{\gamma })$ and the relative optomechanical coupling strength of the jth element to the field. Note that, because of the normalization of the g_j, the total heat flow through the array is independent of the system size⁹ :

$$\begin{eqnarray}&&{{J}_{{\rm m}}}=2\omega \gamma \eta (n-\bar{n})=2\omega \gamma (n-n^{\prime} )\end{eqnarray} \tag{ 45 }$$

where $n^{\prime} =(1-\eta )n+\eta \bar{n}$ is the final occupation number of the collective mode ${{\tilde{b}}_{1}}$. Let us also note that, while the averaged heat flow per element, $\bar{J}\frac{1}{N}{{\sum }_{j}}{{J}_{j}}={{J}_{{\rm m}}}/N$, scales as the inverse of the length of the arrays, as expected from the Fourier law [1], the local heat flow depends on the form of the individual optomechanical coupling. To take an example, for a field whose wavelength is chosen such that the whole array is 'transmissive', the g_j can be shown to have a sinusoidal dependence with the element position in the array [20, 21]

$$\begin{eqnarray}&&{{g}_{j}}=g\sqrt{\frac{2}{N}}{\rm sin} \left( 2\pi \frac{j-1/2}{N} \right)\quad (N\gt 2).\end{eqnarray} \tag{ 46 }$$

One sees from equation (44) that the currents flowing through the different oscillators have quite different behaviours in the large-N limit. Indeed, as shown in figure 3, at the extremities (or the centre) of the array, one has ${{J}_{j}}\sim 1/{{N}^{3}}$ as $N\to \infty$, while for $j\sim N/4$ or $3N/4$, one obtains ${{J}_{j}}\sim 1/N$. Choosing the form of the optomechanical couplings thus offers some freedom in tuning the heat flow through individual elements.

Zoom In Zoom Out Reset image size

4.3.3. Experimental implementation

The generic features discussed previously could in principle be observed in a broad range of optomechanical systems: in the optical domain, these could be arrays of flexible membranes [48], toroidal cavities with indentations [49], optomechanical crystals [50], or ensembles of cold atoms in optical cavities [22, 51]; in the microwave domain, micromechanical elements coupled to superconducting microwave resonator fields [62, 52] are one such possibility. Depending on the system considered, the measurement of the oscillators' thermal occupancy can be performed optically via sideband thermometry [53, 54] or collective mode readout [20, 51], or electrically by functionalizing the elements [55–57] or by coupling the mechanical elements to additional microwave resonators or to artificial atoms, as demonstrated in, e.g., [58, 59]. The two-oscillator scenario could be for instance realized using two micromirrors [60, 61] in the double-cavity geometry of [45] in which the motion of the centre-of-mass and relative-motion modes can be addressed and readout by two optical fields appropriately detuned from the cavity resonance. Alternatively, one could use a pair of partially transmitting, flexible membranes positioned in an optical cavity driven by optical fields with specific wavelengths [20].

The variety of experimental implementations available implies that a large swathe of parameter space is experimentally accessible. For example, in an electromechanical implementation as the one in [62], with $\Omega =2\pi \times 7$ GHz, $\kappa =2\pi \times 6$ MHz, $-\Delta =\omega =2\pi \times 32$ MHz, and γ typically three orders of magnitude smaller than ω [63], one would expect ${{J}_{{\rm m}}}\sim {{10}^{-19}}$ W ($\sim 7\times {{10}^{6}}$ phonons/s) at 25 mK, and assuming efficient optomechanical cooling (${{n}^{\prime }}\ll n$). For optomechanical implementations with silicon nitride membranes, $\Omega \sim 2\pi \times 282$ THz, $\kappa \sim 2\pi \times 1.0$ MHz, $-\Delta =\omega =2\pi \times 136$ kHz, and $\gamma /\omega \sim {{10}^{-6}}$ [64], we find ${{J}_{{\rm m}}}\sim {{10}^{-22}}$ W ($\sim {{10}^{6}}$ phonons/s) at 4 K.

We have investigated the transport of heat in a system of N quantum oscillators coupled to common and independent baths and derived analytical expressions for the steady state occupancy and heat flow. The obtained results are, among others, relevant in the context of optomechanical arrays where by choosing the form of the coupling between the optical field and the mechanics one can engineer effective couplings and baths for the mechanical oscillators, and thereby tune the heat flow through individual elements. While the present work focussed on the situation of an optical field in a coherent state—a common oscillator coupled to a zero-temperature bath—the same approach could be used to tackle the case of optomechanical interactions with (Gaussian) fields exhibiting nonclassical correlations, such as squeezing [45, 46].

This work was supported by the COST Action MP1209 'Thermodynamics in the Quantum Regime', and the Danish Council for Independent Research (Sapere Aude program).

The purpose of this appendix is to introduce to the unfamiliar reader the main steps involved in the adiabatic elimination of a fast system, in our case an optical field described by the annihilation operator $\hat{a}$, from a coupled system involving the optical field and a mechanical mode, described by the annihilation operator $\hat{b}$. Our treatment can easily be generalized to many mechanical modes and follows very closely the exposition in the Supplemental Information of [25]. We refer the interested reader to alternative expositions, e.g., in [47], for further detail. To start, we divide the total Hamiltonian of the system into three components. First is the free Hamiltonian

$$\begin{eqnarray}&&{{\hat{H}}_{{\rm free}}}=-\Delta {{\hat{a}}^{\dagger }}\hat{a}+\omega {{\hat{b}}^{\dagger }}{{\hat{b}}^{{}}},\end{eqnarray} \tag{ 47 }$$

which we write in the rotating from of the driving field (introduced below), such that $\Delta {{\Omega }_{{\rm L}}}-\Omega$ is the difference between the resonance frequency of mode 'a' and its driving field. Next, we must consider the Hamiltonian that couples the two subsystems,

$$\begin{eqnarray}&&{{\hat{H}}_{\operatorname{int}}}=-g\left( \hat{b}+{{{\hat{b}}}^{\dagger }} \right){{\hat{a}}^{\dagger }}\hat{a}.\end{eqnarray} \tag{ 48 }$$

Finally, we add the Hamiltonian that acts to drive the optical field to a non-zero mean value,

$$\begin{eqnarray}&&{{\hat{H}}_{{\rm d}}}={\rm i}\sqrt{2\kappa }{{\alpha }_{{\rm in}}}\left( {{{\hat{a}}}^{\dagger }}-\hat{a} \right);\end{eqnarray} \tag{ 49 }$$

here, the (coherent) driving field amplitude ${{\alpha }_{{\rm in}}}$ is taken to be real for simplicity. The compound system is also acted on by dissipation, which we concisely model using the superoperators

$$\begin{eqnarray}&&D[\hat{c}]\varrho 2\hat{c}\varrho {{\hat{c}}^{\dagger }}-{{\hat{c}}^{\dagger }}\hat{c}\varrho -\varrho {{\hat{c}}^{\dagger }}\hat{c},\end{eqnarray} \tag{ 50 }$$

and

$$\begin{eqnarray}&&{{D}_{{\rm th}}}[\hat{c}](n+1)D[\hat{c}]+nD[{{\hat{c}}^{\dagger }}],\end{eqnarray} \tag{ 51 }$$

where n gives the mean number of excitations present in the bath that $\hat{c}$ is coupled to. We now label the (amplitude) decay rates of $\hat{a}$ and $\hat{b}$ by κ and γ, respectively, such that the full master equation of the system takes the form

$$\begin{eqnarray}&&\dot{\varrho }=-{\rm i}\left[ {{{\hat{H}}}_{{\rm free}}}+{{{\hat{H}}}_{\operatorname{int}}}+{{{\hat{H}}}_{{\rm d}}},\varrho \right]+\kappa D[\hat{a}]\varrho +\gamma {{D}_{{\rm th}}}\left[ {\hat{b}} \right]\varrho .\end{eqnarray} \tag{ 52 }$$

We now proceed by defining two simplifying unitary transformations; we shift both $\hat{a}$ and $\hat{b}$ by an as-yet undetermined complex number:

$$\begin{eqnarray}&&\hat{U}_{{\rm a}}^{\dagger }\hat{a}\hat{U}_{{\rm a}}^{\dagger }=\hat{a}+\alpha \;{\rm and}\;\hat{U}_{{\rm m}}^{\dagger }\hat{b}\hat{U}_{{\rm m}}^{\dagger }=\hat{b}+\beta .\end{eqnarray} \tag{ 53 }$$

We determine α and β by demanding that they act to cancel out the contribution of ${{\hat{H}}_{{\rm d}}}$ in the master equation, that $\hat{a}$ and $\hat{b}$ are shifted to have zero mean value. Furthermore, we suppose that $|\alpha |\gg 1$, such that we can neglect the term in the transformed ${{\hat{H}}_{\operatorname{int}}}$ that is proportional to ${{\hat{a}}^{\dagger }}\hat{a}$. Finally, we obtain a transformed Hamiltonian that is quadratic in the operators:

$$\begin{eqnarray}&&{{\hat{H}}_{{\rm trans}}}-\Delta {{\hat{a}}^{\dagger }}\hat{a}+\omega {{\hat{b}}^{\dagger }}\hat{b}-g\left( \hat{b}+{{{\hat{b}}}^{\dagger }} \right)\left( {{\alpha }^{*}}\hat{a}+\alpha {{{\hat{a}}}^{\dagger }} \right).\end{eqnarray} \tag{ 54 }$$

We have further redefined Δ in order to absorb a shift in the cavity frequency introduced by β. It is now convenient to consider the three terms making up ${{\hat{H}}_{{\rm trans}}}$ separately. In order to do this, we define the system ${{\hat{H}}_{{\rm s}}}\omega {{\hat{b}}^{\dagger }}\hat{b}$, bath ${{\hat{H}}_{{\rm b}}}-\Delta {{\hat{a}}^{\dagger }}\hat{a}$, and interaction ${{\hat{H}}_{{\rm i}}}-g\left( \hat{b}+{{{\hat{b}}}^{\dagger }} \right)\left( {{\alpha }^{*}}\hat{a}+\alpha {{{\hat{a}}}^{\dagger }} \right)$ Hamiltonians. We also define the associated Liouvillians

$$\begin{eqnarray}&&{{\mathcal{L}}_{{\rm s}}}\varrho -{\rm i}\left[ {{{\hat{H}}}_{{\rm s}}},\varrho \right]+\gamma {{D}_{{\rm th}}}\left[ {\hat{b}} \right]\varrho ,\end{eqnarray} \tag{ 55 }$$

$$\begin{eqnarray}&&{{\mathcal{L}}_{{\rm b}}}\varrho -{\rm i}\left[ {{{\hat{H}}}_{{\rm b}}},\varrho \right]+\gamma D[\hat{a}]\varrho ,\end{eqnarray} \tag{ 56 }$$

and

$$\begin{eqnarray}&&{{\mathcal{L}}_{{\rm i}}}\varrho -{\rm i}\left[ {{{\hat{H}}}_{{\rm i}}},\varrho \right],\end{eqnarray} \tag{ 57 }$$

respectively. With this notation, the transformed master equation acquires the simple form $\dot{\varrho }=\left( {{\mathcal{L}}_{{\rm s}}}+{{\mathcal{L}}_{{\rm b}}}+{{\mathcal{L}}_{{\rm i}}} \right)\varrho$. Our goal in this section is to derive an effective equation of motion for the density matrix describing the mechanical subsystem, with mode 'a' eliminated from it. In common with other such eliminations, our results will only be valid to lowest order in the coupling strength g.

The approach we follow uses projection operators to effectively project the master equation and achieve this elimination. The first projection operator we require takes the form

$$\begin{eqnarray}&&\mathbb{P}\varrho {{\rho }_{{\rm a}}}\;\otimes \;{\rm T}{{{\rm r}}_{{\rm a}}}\varrho ,\end{eqnarray} \tag{ 58 }$$

where ${{\rho }_{{\rm a}}}$ is the vacuum state for mode 'a' and corresponds to the solution of the master equation for mode 'a' alone. We define $\mathbb{I}$ as the identity operator and a further projection operator $\mathbb{Q}=\mathbb{I}-\mathbb{P}$. Following [47, section 5.2.1] we note that

$$\begin{eqnarray}&&\mathbb{P}{{\mathcal{L}}_{{\rm i}}}\mathbb{P}=0,{\rm and}\end{eqnarray} \tag{ 59 }$$

$$\begin{eqnarray}&&\mathbb{Q}\left( {{\mathcal{L}}_{{\rm s}}}+{{\mathcal{L}}_{{\rm b}}} \right)=\left( {{\mathcal{L}}_{{\rm s}}}+{{\mathcal{L}}_{{\rm b}}} \right)\mathbb{Q}.\end{eqnarray} \tag{ 60 }$$

These projection operators are used to project the full master equation, resulting in the two differential equations

$$\begin{eqnarray}&&\mathbb{P}\dot{\varrho }=\mathbb{P}\left( {{\mathcal{L}}_{{\rm s}}}+{{\mathcal{L}}_{{\rm b}}} \right)\mathbb{P}\varrho +\mathbb{P}{{\mathcal{L}}_{{\rm i}}}\mathbb{Q}\varrho ,{\rm and}\end{eqnarray} \tag{ 61 }$$

$$\begin{eqnarray}&&\mathbb{Q}\dot{\varrho }=\mathbb{Q}\left( {{\mathcal{L}}_{{\rm s}}}+{{\mathcal{L}}_{{\rm b}}}+{{\mathcal{L}}_{{\rm i}}} \right)\mathbb{Q}\varrho +\mathbb{Q}{{\mathcal{L}}_{{\rm i}}}\mathbb{P}\varrho .\end{eqnarray} \tag{ 62 }$$

We now invoke the weak-coupling approximation, assuming g to be small. We can then formally write

$$\begin{eqnarray}&&\mathbb{Q}\varrho (t)={{{\rm e}}^{\mathbb{Q}({{\mathcal{L}}_{{\rm s}}}+{{\mathcal{L}}_{{\rm b}}}+{{\mathcal{L}}_{{\rm i}]}})(t-{{t}_{0}})}}\mathbb{Q}\varrho ({{t}_{0}})+\int _{{{t}_{0}}}^{t}{\rm d}\tau {{{\rm e}}^{\mathbb{Q}({{\mathcal{L}}_{{\rm s}}}+{{\mathcal{L}}_{{\rm b}}}+{{\mathcal{L}}_{{\rm i}}})(t-\tau )}}\mathbb{Q}{{\mathcal{L}}_{{\rm i}}}\mathbb{P}\varrho (\tau ).\end{eqnarray} \tag{ 63 }$$

Next, we substitute the lowest-order approximation

$$\begin{eqnarray}&&\mathbb{P}\varrho (\tau )={{{\rm e}}^{\mathbb{P}({{\mathcal{L}}_{{\rm s}}}+{{\mathcal{L}}_{{\rm b}}})(\tau -t)}}\mathbb{P}\varrho (t)\end{eqnarray} \tag{ 64 }$$

in the integrand, and we take the initial time ${{t}_{0}}\to -\infty$. Because of our lowest-order-in-g approximation, we ignore the ${{\mathcal{L}}_{{\rm i}}}$ in the exponent, yielding

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \mathbb{P}\dot{\varrho } & \approx & \mathbb{P}\left( {{\mathcal{L}}_{{\rm s}}}+{{\mathcal{L}}_{{\rm b}}} \right)\mathbb{P}\varrho \\ {} & {} & +\int _{0}^{\infty }{\rm d}\tau \mathbb{P}{{\mathcal{L}}_{{\rm i}}}{{{\rm e}}^{\mathbb{Q}({{\mathcal{L}}_{{\rm s}}}+{{\mathcal{L}}_{{\rm b}}})\tau }}\mathbb{Q}{{\mathcal{L}}_{{\rm i}}}{{{\rm e}}^{-\mathbb{P}({{\mathcal{L}}_{{\rm s}}}+{{\mathcal{L}}_{{\rm b}}})\tau }}\mathbb{P}\varrho ({\rm t}). \\ \end{array}\end{eqnarray} \tag{ 65 }$$

Next, we notice that $[\mathbb{Q},{{\mathcal{L}}_{{\rm s}}}+{{\mathcal{L}}_{{\rm b}}}]=0$, so that we may simplify the above differential equation to

$$\begin{eqnarray}&&\mathbb{P}\dot{\varrho }\approx \mathbb{P}\left( {{\mathcal{L}}_{{\rm s}}}+{{\mathcal{L}}_{{\rm b}}} \right)\mathbb{P}\varrho +\int _{0}^{\infty }{\rm d}\tau \mathbb{P}{{\mathcal{L}}_{{\rm i}}}\mathbb{Q}{{{\rm e}}^{({{\mathcal{L}}_{{\rm s}}}+{{\mathcal{L}}_{{\rm b}}})\tau }}{{\mathcal{L}}_{{\rm i}}}{{{\rm e}}^{-\mathbb{P}({{\mathcal{L}}_{{\rm s}}}+{{\mathcal{L}}_{{\rm b}}})\tau }}\left( {{\rho }_{{\rm a}}}\;\otimes \;{\rm T}{{{\rm r}}_{{\rm a}}}\varrho \right)\end{eqnarray} \tag{ 66 }$$

$$\begin{eqnarray}&&=\mathbb{P}\left( {{\mathcal{L}}_{{\rm s}}}+{{\mathcal{L}}_{{\rm b}}} \right)\mathbb{P}\varrho +\int _{0}^{\infty }{\rm d}\tau \mathbb{P}{{\mathcal{L}}_{{\rm i}}}\mathbb{Q}{{{\rm e}}^{({{\mathcal{L}}_{{\rm s}}}+{{\mathcal{L}}_{{\rm b}}})\tau }}{{\mathcal{L}}_{{\rm i}}}\left( {{\rho }_{{\rm a}}}\;\otimes \;{{{\rm e}}^{-{{\mathcal{L}}_{{\rm s}}}\tau }}{\rm T}{{{\rm r}}_{{\rm a}}}\varrho \right).\end{eqnarray} \tag{ 67 }$$

We can now trace out mode 'a', thereby obtaining an approximate master equation for the reduced density matrix $\rho {\rm T}{{{\rm r}}_{{\rm a}}}\varrho$,

$$\begin{eqnarray}&&\dot{\rho }={{\mathcal{L}}_{{\rm s}}}{{\rho }_{{\rm b}}}+{\rm T}{{{\rm r}}_{{\rm a}}}\int _{0}^{\infty }{\rm d}\tau {{\mathcal{L}}_{{\rm i}}}\mathbb{Q}{{{\rm e}}^{({{\mathcal{L}}_{{\rm s}}}+{{\mathcal{L}}_{{\rm b}}})\tau }}{{\mathcal{L}}_{{\rm i}}}\left( {{\rho }_{{\rm a}}}\;\otimes \;{{{\rm e}}^{-{{\mathcal{L}}_{{\rm s}}}\tau }}\rho \right).\end{eqnarray} \tag{ 68 }$$

The integral can be evaluated by using the definition of ${{\mathcal{L}}_{{\rm i}}}$ given above. Upon applying the rotating-wave approximation, where quickly-rotating terms are neglected, and performing some algebra, the standard form

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {\dot{\rho }} & = & {{\mathcal{L}}_{{\rm s}}}\rho -{\rm i}{{g}^{2}}|\alpha {{|}^{2}}{\rm Im}S(\omega )+S(-\omega )\left[ {{{\hat{b}}}^{\dagger }}\hat{b},\rho \right] \\ {} & {} & +{{g}^{2}}|\alpha {{|}^{2}}\left\{ \operatorname{Re}S(\omega )D\left[ {\hat{b}} \right]\rho +\operatorname{Re}S(-\omega )D\left[ {{{\hat{b}}}^{\dagger }} \right]\rho \right\} \\ \end{array}\end{eqnarray} \tag{ 69 }$$

is recovered. In this expression, we have defined the spectral density of the cavity field to zeroth order in g, which can be expressed as

$$\begin{eqnarray}&&S(\omega )\int _{0}^{\infty }{\rm d}\tau {{{\rm e}}^{{\rm i}\omega \tau }}\langle \hat{a}(t+\tau ){{\hat{a}}^{\dagger }}(t)\rangle =-\frac{1}{{\rm i}(\Delta +\omega )-\kappa }.\end{eqnarray} \tag{ 70 }$$

By comparing equation (69) to the standard master equation for a damped harmonic oscillator, we can identify an effective Hamiltonian

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{{\hat{H}}}_{{\rm eff}}} & = & {{{\hat{H}}}_{{\rm s}}}+{{g}^{2}}|\alpha {{|}^{2}}{\rm Im}S(\omega )+S(-\omega ){{{\hat{b}}}^{\dagger }}\hat{b} \\ {} & = & {{{\hat{H}}}_{{\rm s}}}+\frac{2{{g}^{2}}|\alpha {{|}^{2}}\Delta \left( {{\Delta }^{2}}-{{\omega }^{2}}+{{\kappa }^{2}} \right)}{{{\left( {{\Delta }^{2}}-{{\omega }^{2}}-{{\kappa }^{2}} \right)}^{2}}+{{(2\Delta \kappa )}^{2}}}{{{\hat{b}}}^{\dagger }}\hat{b}, \\ \end{array}\end{eqnarray} \tag{ 71 }$$

as well as the effective cooling and heating Liouvillians

$$\begin{eqnarray}&&{{\mathcal{L}}_{{\rm cool}}}\left[ (n+1)\gamma +\frac{{{g}^{2}}|\alpha {{|}^{2}}\kappa }{{{(\Delta +\omega )}^{2}}+{{\kappa }^{2}}} \right]D\left[ {\hat{b}} \right],\end{eqnarray} \tag{ 72 }$$

and

$$\begin{eqnarray}&&{{\mathcal{L}}_{{\rm heat}}}\left[ n\gamma +\frac{{{g}^{2}}|\alpha {{|}^{2}}\kappa }{{{(\Delta -\omega )}^{2}}+{{\kappa }^{2}}} \right]D\left[ {{{\hat{b}}}^{\dagger }} \right].\end{eqnarray} \tag{ 73 }$$

The equations of motion for the second-order moments obtained from equation (37) are

$$\begin{eqnarray*}\begin{array}{ccccccccccccccc} {{\partial }_{t}}\langle \hat{x}_{1}^{2}\rangle & = & \omega \langle \{{{{\hat{x}}}_{1}},{{{\hat{p}}}_{1}}\}\rangle , \\ {{\partial }_{t}}\langle \hat{x}_{2}^{2}\rangle & = & \omega \langle \{{{{\hat{x}}}_{2}},{{{\hat{p}}}_{2}}\}\rangle , \\ {{\partial }_{t}}\langle {{{\hat{x}}}_{1}}{{{\hat{x}}}_{2}}\rangle & = & {{\partial }_{t}}\langle {{{\hat{x}}}_{2}}{{{\hat{x}}}_{1}}\rangle =\omega (\langle {{{\hat{x}}}_{1}}{{{\hat{p}}}_{2}}\rangle +\langle {{{\hat{x}}}_{2}}{{{\hat{p}}}_{1}}\rangle ), \\ {{\partial }_{t}}\langle \{{{{\hat{x}}}_{1}},{{{\hat{p}}}_{1}}\}\rangle & = & 2\omega (\langle \hat{p}_{1}^{2}\rangle -\langle \hat{x}_{1}^{2}\rangle )-2\Lambda \langle {{{\hat{x}}}_{1}}{{{\hat{x}}}_{2}}\rangle -(2\gamma +\bar{\gamma })\langle \{{{{\hat{x}}}_{1}},{{{\hat{p}}}_{1}}\}\rangle \\ {} & {} & -2\bar{\gamma }\langle {{{\hat{x}}}_{1}}{{{\hat{p}}}_{2}}\rangle , \\ {{\partial }_{t}}\langle \{{{{\hat{x}}}_{2}},{{{\hat{p}}}_{2}}\}\rangle & = & 2\omega (\langle \hat{p}_{2}^{2}\rangle -\langle \hat{x}_{2}^{2}\rangle )-2\Lambda \langle {{{\hat{x}}}_{1}}{{{\hat{x}}}_{2}}\rangle -(2\gamma +\bar{\gamma })\langle \{{{{\hat{x}}}_{2}},{{{\hat{p}}}_{2}}\}\rangle \\ {} & {} & -2\bar{\gamma }\langle {{{\hat{p}}}_{1}}{{{\hat{x}}}_{2}}\rangle , \\ {{\partial }_{t}}\langle \{{{{\hat{x}}}_{1}},{{{\hat{p}}}_{2}}\}\rangle & = & 2\omega (\langle {{{\hat{p}}}_{1}}{{{\hat{p}}}_{2}}\rangle -\langle {{{\hat{x}}}_{1}}{{{\hat{x}}}_{2}}\rangle )-2\Lambda \langle \hat{x}_{1}^{2}\rangle -2(2\gamma +\bar{\gamma })\langle {{{\hat{x}}}_{1}}{{{\hat{p}}}_{2}}\rangle \\ {} & {} & -\bar{\gamma }\langle \{{{{\hat{x}}}_{1}},{{{\hat{p}}}_{1}}\}\rangle , \\ {{\partial }_{t}}\langle \{{{{\hat{x}}}_{2}},{{{\hat{p}}}_{1}}\}\rangle & = & 2\omega (\langle {{{\hat{p}}}_{1}}{{{\hat{p}}}_{2}}\rangle -\langle {{{\hat{x}}}_{1}}{{{\hat{x}}}_{2}}\rangle )-2\Lambda \langle \hat{x}_{2}^{2}\rangle -2(2\gamma +\bar{\gamma })\langle {{{\hat{x}}}_{2}}{{{\hat{p}}}_{1}}\rangle \\ {} & {} & -\bar{\gamma }\langle \{{{{\hat{x}}}_{2}},{{{\hat{p}}}_{2}}\}\rangle , \\ {{\partial }_{t}}\langle \hat{p}_{1}^{2}\rangle & = & 7-\omega \langle \{{{{\hat{x}}}_{1}},{{{\hat{p}}}_{1}}\}\rangle -2\Lambda \langle {{{\hat{x}}}_{2}}{{{\hat{p}}}_{1}}\rangle -2(2\gamma +\bar{\gamma })\langle \hat{p}_{1}^{2}\rangle -2\bar{\gamma }\langle {{{\hat{p}}}_{1}}{{{\hat{p}}}_{2}}\rangle \\ {} & {} & +\left[ 2\gamma (2{{n}_{1}}+1)+\bar{\gamma }(2\bar{n}+1) \right], \\ {{\partial }_{t}}\langle \hat{p}_{2}^{2}\rangle & = & -\omega \langle \{{{{\hat{x}}}_{2}},{{{\hat{p}}}_{2}}\}\rangle -2\Lambda \langle {{{\hat{x}}}_{1}}{{{\hat{p}}}_{2}}\rangle -2(2\gamma +\bar{\gamma })\langle \hat{p}_{2}^{2}\rangle -2\bar{\gamma }\langle {{{\hat{p}}}_{1}}{{{\hat{p}}}_{2}}\rangle \\ {} & {} & +\left[ 2\gamma (2{{n}_{2}}+1)+\bar{\gamma }(2\bar{n}+1) \right], \\ \end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}\begin{array}{ccccccccccccccc} {{\partial }_{t}}\langle {{{\hat{p}}}_{1}}{{{\hat{p}}}_{2}}\rangle & = & {{\partial }_{t}}\langle {{{\hat{p}}}_{2}}{{{\hat{p}}}_{1}}\rangle =-\omega (\langle {{{\hat{x}}}_{1}}{{{\hat{p}}}_{2}}\rangle +\langle {{{\hat{p}}}_{1}}{{{\hat{x}}}_{2}}\rangle )-\Lambda (\{{{{\hat{x}}}_{1}}{{{\hat{p}}}_{1}}\}+\{{{{\hat{p}}}_{2}}{{{\hat{x}}}_{2}}\}) \\ {} & {} & -2(2\gamma +\bar{\gamma })\langle {{{\hat{p}}}_{1}}{{{\hat{p}}}_{2}}\rangle -\bar{\gamma }(\langle \hat{p}_{1}^{2}\rangle +\langle \hat{p}_{2}^{2}\rangle )+\bar{\gamma }(2\bar{n}+1). \\ \end{array}\end{eqnarray*}$$

Recall also that $\langle [{{\hat{x}}_{j}},{{\hat{p}}_{k}}]\rangle ={\rm i}{{\delta }_{j,k}}$, which gives ${{\partial }_{t}}\langle [{{\hat{x}}_{j}},{{\hat{p}}_{k}}]\rangle =0$ $(j,k=1,2)$.