Heat transport in harmonic oscillator systems with correlated baths: Application to optomechanical arrays

We investigate the transport of phonons between $N$ harmonic oscillators in contact with independent thermal baths and coupled to a common oscillator, and derive an expression for the steady state heat flow between the oscillators in the weak coupling limit. We apply these results to an optomechanical array consisting of a pair of mechanical resonators coupled to a single quantised electromagnetic field mode by radiation pressure as well as to thermal baths with different temperatures. In the weak coupling limit this system is shown to be equivalent to two mutually-coupled harmonic oscillators in contact with an effective common thermal bath in addition to their independent baths. The steady state occupation numbers and heat flows are derived and discussed in various regimes of interest.


Introduction
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][2][3][4][5][6][7][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-ofequilibrium 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][12][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][25][26][27][28][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][32][33][34][35][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.
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-ofequilibrium system. Finally, we conclude in section 5 by surveying our results and putting them in the context of possible future work.

2.
Heat flow for N oscillators coupled to a common oscillator and independent thermal baths 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 Figure 1. (a) N oscillators, each in contact with its own independent thermal bath and coupled to a common oscillator, which is in contact with its own bath. (b) Equivalent model system for an array with N = 2: after adiabatic elimination of the common oscillator, the oscillators are effectively mutually coupled and in contact with both their initial thermal bath and a common bath. density matrix for the + N ( 1)-partite system. The master equation for ϱ can be written ( =  1) : where Ĥ is the Hamiltonian governing the evolution, and bˆj the annihilation operator of the jth oscillator ( ⩽ ⩽ j N 1 ), γ j the coupling rate of the jth oscillator to its bath, whose mean occupation number is n j , and with â the annihilation operator of the final harmonic oscillator, and κ its coupling constant to its own bath, assumed at zero-temperature. We assume that Ĥ has the following form [25]: 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, â can be adiabatically eliminated from the dynamics 4 . 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 5 is given by [5]  ϱ = ( ) We now make use of our adiabatic elimination procedure to write, in steady state and to lowest order in the coupling constants g j , ϱ ρ ρ = ⊗ ss a , where ρ a is the steady-state density matrix for the + N ( 1)stharmonic oscillator. Let us first take the trace with respect to this mode 6 : Tr TrˆˆˆˆˆˆˆT rˆˆˆˆˆˆˆˆ.
Following the usual methodology, and under the assumption that the oscillator described by â is driven into a coherent state, the state described by ρ a has been shifted to describe a zero-mean vacuum state, so that . 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 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 4 The assumptions of the adiabatic elimination procedure are mainly (i) that one subsystem evolves on a much faster time-scale than the rest of the system, (ii) that the coupling between the two is weak, and (iii) that the total density matrix is approximately a tensor product between the fast and slow subsystems. 5 We use the convention where positive heat flow corresponds to heat flowing into the system. heat flow: 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.

Case of identical baths
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 γ γ = j and occupation number = n n j . To continue, it is convenient to introduce a normal-mode basis. Define ( ) forms an orthonormal basis 7 . We use this 'collective' basis to define a new set of N normalized harmonic oscillator annihilation operators b j ( ⩽ ⩽ j N 1 ): where each  j is defined analogously to  j but with bˆj replaced by b j . The Hamiltonian can be expressed as the sum of Hamiltonians for − N 1 uncoupled and two linearly-coupled oscillators: Note that it is only because the thermal baths all have the same temperature that the Liouvillians L j are diagonal in the new basis. As shown in appendix A, the adiabatic elimination of â and application of the rotating-wave approximation leads to a shift of the oscillator 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 The reduced master equation (13) has a steady-state solution ρ ss given by a tensor product of N thermal states ρ j .
We define ′ = 〈 〉 n b b1 † 1 to be the occupation number of the state for j = 1, and = 〈 〉 n b bj j † for ⩾ j 2. Writing the heat flow in the normal-mode basis yields 7 Since we define g j such that they are all real, we can assume that the entire basis set is composed of entirely real vectors.
Because of the uncorrelated nature of the steady state, we can immediately write that where we have also taken the sum out of the parentheses by exploiting the properties of the orthonormal matrix G. By the very nature of the steady state, however, Therefore, the heat flowing into or out of the mechanical subsystem is simply which is nonzero because we have traced out the + N ( 1)stoscillator. Let us note that − ′ → n n 0 when → g 0. The heat flowing into or out of this oscillator must therefore be 4. Application to optomechanical arrays 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 1 and n 2 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 8 . Under these conditions the Hamiltonian and Liouvillian of the system are given by We have defined the cavity field detuning Δ Ω where Ω L is the frequency of the driving field, assumed to be monochromatic, and have transformed our system to a frame rotating at the frequency Ω L . The enhanced optomechanical coupling rate g and the Liouvillians  1 (2) and  a are defined as in equations (2) and (4), respectively. Introducing the relative and centre-of-mass modes these can be recast as r c r , c a with  r,c a Liouvillian including the correlations between the rotated baths. 8 Note that the opposite situation can be realized in, e.g., a 'transmissive' configuration [20] or in the double-cavity geometry of [45].
Adiabatically eliminating the cavity field in the weak optomechanical coupling regime leads to an effective mechanical frequency ω ω Λ ′ = + , due to the optical spring effect, as well as an effective damping rate γ γ γ ′ = +¯, for the centre-of-mass mode (see appendix A), where The field fluctuations also give rise to an additional effective bath, coupling to the centre-of-mass mode only, whose Liouvillian has the form Assuming a red-detuned cavity field (Δ < 0) for which β β > + − , this term describes the coupling with a thermal bath with coupling rate and occupation number defined by The effective Hamiltonian and Liouvillian after adiabatic elimination then read eff r c r,c or, going back to the bare mechanical basis, eff 1 2 In the regime of interest, Λ ω ≪ , so that we may approximate 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  1,2 and  , 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.

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]: 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 = 1, 2). We assume the baths to be Markovian, the independent bath temperatures being given by n j , and the common bath temperature by n. In order to obtain an effective master equation in the limit γ γ Λ ω ≪ ≪ ,¯, we employ these assumptions to follow the procedure of [4], which we will not detail explicitly but finally yields Under the rotating-wave approximation, where terms of the form b b1 2 and b b1 † 2 † 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, and thereby calculate the heat flow in this system from (9). One gets ( 1) 2¯2¯( 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.

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, κ ω ≫ . In this regime, and for Δ κ ∼ − , equations (27) and (28) give In contrast, the regime of large coupling to the effective common bath can be achieved in the resolved sideband regime of optomechanics, ω κ (i) In the large coupling regime Λ γ γ ≫ ,¯, one has which we can understand in two different limits. (a) When the coupling to the independent baths is larger than that to the common bath (γ γ ≫¯), the mean independent bath temperature + n n ( )2 1 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 j j This is consistent with the limit of radiation pressure cooling in optomechanics [41][42][43][44], which predicts that the (centre-of-mass) mode coupled to the field is cooled down to n, γ γ ≫ . Since the uncoupled (relative motion) mode's occupancy remains n, this means that, in the bare basis, ′ → + n n n (¯) 2 j when γ γ ≫ .
As an illustration, figure 2 shows the individual heat flows J 1 and J 2 , as well as the total heat flow through the mechanics + J J 1 2 , as a function of Λ γ and γ γ , 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 γ, 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.

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 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 η γ γ γ = + (¯) 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 9 : , 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] One sees from equation (44)

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][56][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 32MHz, and γ typically three orders of magnitude smaller than ω [63]  The averaged heat flow per element J (mustard) is also shown. 9 We remark here that any intrinsic dependence of g itself on the system size can be counteracted by varying the driving strength appropriately.

Conclusion
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 zerotemperature 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]. . 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 These projection operators are used to project the full master equation, resulting in the two differential equations We now invoke the weak-coupling approximation, assuming g to be small. We can then formally write Next, we substitute the lowest-order approximation in the integrand, and we take the initial time → −∞ t 0 . Because of our lowest-order-in-g approximation, we ignore the  i in the exponent, yielding