Cavity-mediated collective laser-cooling of a non-interacting atomic gas inside an asymmetric trap to very low temperatures

Abstract

In this paper, we identify a many-particle phonon expectation value $\zeta$ with the ability to induce collective dynamics in a non-interacting atomic gas inside an optical cavity. We then propose to utilise this expectation value to enhance the laser cooling of many atoms through a cyclic two-stage process which consists of cooling and displacement stages. During cooling stages, short laser pulses are applied. These use $\zeta$ as a resource and decrease the vibrational energy of the atomic gas by a fixed amount. Subsequent displacement stages use the asymmetry of the trapping potential to replenish the many-particle phonon expectation value $\zeta$. Alternating both stages of the cooling process is shown to transfer the atomic gas to a final temperature which vanishes in the infinitely-many particle limit.

Keywords:

• Cavity QED
• collective effects
• laser cooling

1. Introduction

Laser sideband cooling allows for cooling single, strongly-confined atomic particles to very low temperatures [1,2]. Its discovery opened the way for experiments which test the foundations of the quantum physics of single particles and have applications ranging from quantum metrology to quantum computing [3]. However, despite imminent quantum technological applications of ultra-cold many-body systems, applying laser cooling to large ensembles of particles or more complex species, like molecules, is not straightforward. For example, the most efficient technique to cool atomic gases to the low temperatures required for Bose Einstein condensation is evaporative cooling, which involves the removal of a significant percentage of particles from the trap. Alternative cooling techniques with the potential of maintaining the initial number of trapped particles therefore receive a lot of attention in the literature [4].

For example, it is hoped that cavity-mediated laser-cooling will enable us to cool large samples of atomic particles collectively and non-destructively to very low temperatures. First indications that this might indeed be the case were found in Paris already more than twenty years ago [5]. More systematic experimental studies of cavity-mediated laser cooling have subsequently been reported by several groups [6–14]. All of these experiments have their respective merits. However, all of them are essentially equivalent to laser-cooling without optical cavities and do not provide a much more efficient cooling mechanism for many particles.

The theory of cavity-mediated laser cooling of free particles was first discussed by Mossberg et al. [15] and Zaugg et al. [16]. Later, Ritsch and collaborators [17,18] and others [19–21] developed semiclassical theories to model cavity-mediated cooling processes. In 1993, Cirac et al. [22] introduced a master equation approach to analyse cavity-mediated laser-cooling in more detail. Since the precision of calculations, which are based on master equations, is easier to control than the precision of semiclassical approximations, this approach has subsequently been used by many authors to emphasize a close analogy between laser-sideband and cavity-mediated laser-cooling of atomic particles [23–28].

Recently, Mishina showed that the cooling dynamics of particles can be affected by trap inhomogeneities [29]. In this paper, we adopt an alternative but conceptually related approach to cool a non-interacting atomic gas inside an optical cavity. Instead of considering a symmetric harmonic trap, we place the particles into an asymmetric trapping potential. The cooling process that we propose here is qualitatively different from other laser cooling techniques [22–28] for which the final energy of the atoms is usually independent of the number N of particles inside the trap. In contrast to this, we predict a final mean phonon number which vanishes in the infinitely-many particle limit. Hence the proposed cooling mechanismhas the potential to transfer an atomic gas to very low temperatures. Notice that the paper by Mishina [29] predicts a similarly strong but less useful dependence of the final temperature of the atoms on the particle number.

Figure 1. Schematic view of the dynamics of the mean phonon number m and the mean phonon expectation value $\zeta$, which is a measure for the average distance of the atoms from the centre of the trapping potential. The proposed cooling process consists of cooling stages interspersed with displacement stages. During each displacement stage, the atoms move freely inside an asymmetric trap and accumulate some distance away from the centre of the trap, thereby accumulating a non-zero value for $\zeta$, while m remains essentially constant. During each cooling stage, a cooling laser with Rabi frequency $\mathrm{\Omega }$ is turned on which returns the atoms to the centre of the trap, thereby returning $\zeta$ to zero while reducing m with a collectively-enhanced cooling rate.

To illustrate our main ideas we first consider an idealised version of the proposed collective cooling mechanism. Our initial toy model is relatively simple and cannot be implemented physically but shows clearly what we aim to achieve, at least approximately, with the experimental setup which we analyse later on in this paper. Suppose N non-interacting atomic particles are so strongly confined that their motion becomes quantised. Moreover we assume in the following that the atomic gas is placed inside an optical cavity and that its system Hamiltonian equals(1) $$\begin{array}{cc}\hfill H_{\mathrm{ideal}}& =\sum _{i=1}^NħJ\left[b_i{c}^{†}+b_i^{†}c\right]\hfill \end{array}$$(1)

with respect to an accordingly chosen interaction picture. Here $b_i$ with $[b_i,b_j^{†}]={\delta }_{i,j}$ and c with $[c,c^{†}]=1$ are the bosonic annihilation operators of a single phonon in the motion of particle i and of a photon inside the resonator, respectively, and J denotes a phonon-photon coupling rate. Moreover, we assume that the field inside the cavity has a finite spontaneous decay rate $\kappa$. Consequently, the density matrix ${\rho }_{\mathrm{I}}$ of phonons and cavity photons evolves in the corresponding interaction picture according to a Markovian master equation in Lindblad form,(2) $$\begin{array}{c}\hfill {\dot{\rho }}_{\mathrm{I}}=-\frac{\mathrm{i}}{ħ}\left[H_{\mathrm{ideal}},{\rho }_{\mathrm{I}}\right]+\frac{\kappa }{2}\left(2c{\rho }_{\mathrm{I}}{c}^{†}-c^{†}c{\rho }_{\mathrm{I}}-{\rho }_{\mathrm{I}}{c}^{†}c\right).\end{array}$$(2)

This master equation can be used to predict the dynamics of expectation values. In the following, we are especially interested in the two expectation values(3) $$\begin{array}{ccc}\hfill m& \equiv \hfill & \hfill \frac{1}{N}\sum _{i=1}^N⟨b_i^{†}{b}_i⟩,\\ \hfill \zeta & \equiv \hfill & \hfill \frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}⟨b_i^{†}{b}_j⟩.\end{array}$$(3)

Both variables are normalised such that they are independent of N in the infinitely many-particle limit. Proceeding as described in Appendix A, one can then show that(4) $$\begin{array}{cc}\hfill \dot{m}=\dot{\zeta }& =-\frac{4N{J}^2}{\kappa }\zeta \hfill \end{array}$$(4)

in case of a relatively large cavity decay rate $\kappa$ and in case of a sufficiently large N,(5) $$\begin{array}{c}\hfill \kappa \gg J\mathrm{and}N\gg 1.\end{array}$$(5)

This implies cooling, if m is at least an approximate measure for the vibrational energy of the trapped particles. The derivation of the above effective rate equations only involves standard approximations which are well justified in the parameter regime given in Equation (5(5) $$\begin{array}{c}\hfill \kappa \gg J\mathrm{and}N\gg 1.\end{array}$$(5) ). In this regime, the term on the right hand side of Equation (4(4) $$\begin{array}{cc}\hfill \dot{m}=\dot{\zeta }& =-\frac{4N{J}^2}{\kappa }\zeta \hfill \end{array}$$(4) ) dominates the time evolution of m and $\zeta$ as long as $\zeta$ is non-zero and sufficiently large.

Notice also that the above-described collective dynamics is the result of the presence of many particles. It is qualitatively different from the dynamics of the individual system components. For example, for $N=1$, Equation (3(3) $$\begin{array}{ccc}\hfill m& \equiv \hfill & \hfill \frac{1}{N}\sum _{i=1}^N⟨b_i^{†}{b}_i⟩,\\ \hfill \zeta & \equiv \hfill & \hfill \frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}⟨b_i^{†}{b}_j⟩.\end{array}$$(3) ) automatically implies $\zeta \equiv 0$. Hence, for a single particle, we always have $\dot{m}=\dot{\zeta }=0$ on the fast time scale. However, for $N\gg 1$ and when $\zeta$ is sufficiently large, the above-described collective dynamics dominates all other dynamics and the mean phonon number m changes relatively rapidly. Suppose $\zeta (0)$ is positive at $t=0$, when the ideal Hamiltonian $H_{\mathrm{ideal}}$ is switched on. Solving Equation (4(4) $$\begin{array}{cc}\hfill \dot{m}=\dot{\zeta }& =-\frac{4N{J}^2}{\kappa }\zeta \hfill \end{array}$$(4) ), we then find that(6) $$\begin{array}{cc}\hfill m(t)& =m(0)-\zeta (0),\hfill \\ \hfill \zeta (t)& =0\hfill \end{array}$$(6)

after a relatively short time $t\gg \kappa /4N{J}^2$. Using this dynamics, it is in principle possible to remove all vibrational energy from the atomic particles and to cool an atomic gas to very low temperatures. To achieve this, the quantum resource $\zeta$, which gets used up in the cooling process, needs to be continuously replenished. This can be done through a cyclic two-step process which consists of cooling and displacement stages, as illustrated in Figure 1. Thecooling stages employ the variable $\zeta$ to rapidly decrease the vibrational energy of the atomic gas. The purpose of the displacement stages is to restore a non-zero and positive phonon expectation value $\zeta$.

Next, let us discuss in more detail, how to replenish $\zeta$, whenever this resource has been used up, without increasing the vibrational energy of the trapped particles. As we shall see below, $\zeta$ is essentially a measure for the mean distance of the atoms from the centre of the trap. One way of solving this task is to employ an asymmetric trapping potential and to let the atoms evolve freely between cooling stages. This is exactly what we propose here. During displacement stages, the atoms oscillate freely within the asymmetric trap and naturally accumulate some distance away from its centre. There is nothing quantum about their motion. The hotter the atoms are, the further they move, while their vibrational energy remains the same.

Figure 2. Schematic view of the proposed experimental setup. This consists of an atomic gas which is externally confined such that the centre of the trap coincides with the node of an optical cavity. The motion should become quantised in the direction of the cavity axis. The trapping potential is assumed to be only approximately harmonic. It should exhibit an asymmetry, which places freely-moving atoms on average slightly away from the centre of the trap. To cool the atomic gas, a series of short laser pulses is applied, which enter the cavity from the side.

Moreover, we need to find a way to realise the system Hamiltonian $H_{\mathrm{ideal}}$ in Equation (1(1) $$\begin{array}{cc}\hfill H_{\mathrm{ideal}}& =\sum _{i=1}^NħJ\left[b_i{c}^{†}+b_i^{†}c\right]\hfill \end{array}$$(1) ), at least to a very good approximation. To do so, we proceed as usual in cavity-mediated laser cooling [23–28]. During each cooling stage, a cooling laser with Rabi frequency $\mathrm{\Omega }$ is applied. The cooling laser should enter the resonator from the side, as illustrated in Figure 2. Moreover, the laser frequency ${\omega }_{\mathrm{L}}$ should coincide with the atomic transition frequency ${\omega }_0$ but should be below the cavity photon frequency ${\omega }_{\mathrm{cav}}$. As we shall see below, best results are obtained when the cavity-atom detuning equal the phonon frequency $\nu$ of the trapped particles,(7) $$\begin{array}{cc}\hfill {\omega }_{\mathrm{cav}}-{\omega }_0& =\nu ,\hfill \end{array}$$(7)

as illustrated in Figure 3. Given this condition, the creation of a cavity photon is most likely accompanied by the loss of a phonon and vice versa, while the atoms remain essentially in a stationary state. Important for the proposed cooling scheme to work is that the atoms are strongly confined. They should be trapped well within an optical wavelength (see e.g. Refs. [30–32]) and the centre of their asymmetric trap should coincide with a node of the cavity field. Moreover, the trap asymmetry should be relatively weak such that the trapping potential is approximately harmonic.

Figure 3. Schematic view of the relevant frequencies of the experimental setup in Figure 2. To maximise the cooling process, the energy of the cavity photons $ħ{\omega }_{\mathrm{cav}}$ should equal the sum of the energy of a single phonon $ħ\nu$ and the energy of an atomic excitation $ħ{\omega }_0$. In addition we assume resonant laser driving of the atoms with frequency ${\omega }_{\mathrm{L}}$ such that ${\omega }_{\mathrm{L}}={\omega }_0$.

For most cavity-mediated laser cooling schemes, it is the ratio between cooling and heating rates which determines the final temperature of the atomic gas (see e.g. Ref. [28]). Any attempt to collectively increase cooling rates usually also results in a collective enhancement of heating rates. As a result, the cooling process might become faster but the final vibrational energy per particle usually remains the same. In contrast to this, for the cooling scheme which we propose here, the time derivative of the expectation value m is always negative. Only cooling processes are collectively enhanced. Heating processes become secondary and the final temperature of the particles is no longer restricted by the system parameters. Alternating displacement and cooling stages, as illustrated in Figure 1, should therefore allows us to transfer the atomic gas to very low temperatures. Cooling only stops when the particles are so cold that they no longer change their position during displacement stages. As we shall see below, the final mean phonon number $m_{\mathrm{final}}$ per particle scales as $1/N^{1/2}$ and tends to zero in the infinitely many-particle limit.

There are five sections in this paper. Section 2 describes the proposed experimental setup and introduces its theoretical model. In Section 3, we analyse the dynamics of the proposed cyclic two-stage cooling process in detail. The asymmetry of the trap and other effects, like detunings and additional interaction terms, which are not present in the above simplified toy-model are explicitly taken into account. We derive effective rate equations for the time evolution of the mean phonon expectation values m and $\zeta$ during cooling and displacement stages. In Section 4, we estimate the general scaling of the final mean phonon number $m_{\mathrm{final}}$ of the proposed cooling scheme. Finally we summarise our findings in Section 5. Some calculations are confined to appendices.

2. Experimental setup and theoretical models

In this section, we discuss the experimental setup in Figure 2 in more detail and introduce the theoretical models which allows for analysing the dynamics of the atomic gas throughout the proposed cooling process.

2.1. Displacement stage

As mentioned already above, the cooling laser is turned off during each displacement stage. When this happens, there is no interaction between vibrational and other degrees of freedom of the atoms. Only the Hamiltonian describing the motion of the particles has to be taken into account when analysing the dynamics of phonon expectation values during this stage. In the following, we assume that the external trap confines the atoms in the direction of the cavity axis. In this case, the Hamiltonian for the vibrational energy of the atomic gas can be written as [28](8) $$\begin{array}{cc}\hfill H_{\mathrm{vib}}& =\sum _{i=1}^N\frac{p_i^2}{2M}+\frac{M{\nu }^2{x}_i^2}{2}+V(x_i),\hfill \end{array}$$(8)

where M is the mass of a single atom. Moreover, $p_i$ is the momentum of particle i in the direction of the cavity axis, $x_i$ denotes its distance from the trap centre and $V(x_i)$ describes corrections to the nearest harmonic trapping potential with phonon frequency $\nu$. To simplify this Hamiltonian, we now introduce bosonic phonon annihilation and creation operators $b_i$ and $b_i^{†}$,(9) $$\begin{array}{ccc}\hfill b_i& \equiv \hfill & \hfill \sqrt{\frac{M\nu }{2ħ}}{x}_i+\mathrm{i}\sqrt{\frac{1}{2Mħ\nu }}{p}_i,\\ \hfill b_i^{†}& \equiv \hfill & \hfill \sqrt{\frac{M\nu }{2ħ}}{x}_i-\mathrm{i}\sqrt{\frac{1}{2Mħ\nu }}{p}_i,\end{array}$$(9)

of particle i with $[b_i,b_i^{†}]=1$. When substituting these operators into Equation (8(8) $$\begin{array}{cc}\hfill H_{\mathrm{vib}}& =\sum _{i=1}^N\frac{p_i^2}{2M}+\frac{M{\nu }^2{x}_i^2}{2}+V(x_i),\hfill \end{array}$$(8) ), neglecting a constant term with no physical consequences and assuming, as an example, a cubic correction $V(x_i)$, the above Hamiltonian $H_{\mathrm{vib}}$ simplifies to(10) $$\begin{array}{cc}\hfill H_{\mathrm{vib}}& =\sum _{i=1}^Nħ\nu b_i^{†}{b}_i+ħ\mu {(b_i+b_i^{†})}^3,\hfill \end{array}$$(10)

where $\mu$ is a constant.

In order to enable us to analytically predict the dynamics of the mean phonon number m and the phonon expectation value $\zeta$, we need to apply certain approximations. Firstly, we assume in the following that the trap is only slightly anharmonic, which implies $\mu \ll \nu$. Moreover, we assume that the trap confines the atoms within an area that is much smaller than the wavelength of the optical resonator. Achieving sub wavelength confinement of atomic particles inside an optical cavity is feasible with current technology [30–32].

2.2. Cooling stage

As illustrated in Figure 1, each displacement stage is followed by a short resonant laser pulse with a (real) Rabi frequency $\mathrm{\Omega }$. The experimental setup, which we need to consider when analysing the cooling stage, is shown in Figure 2. The laser pulse changes the population of the electronic states of the trapped atoms, while affecting their vibrational degrees of freedom and placing photon excitations into the cavity field. To take this into account, we denote the ground state and the excited state of particle i at position $r_i$ by ${|0⟩}_i$ and ${|1⟩}_i$, respectively, while ${\sigma }_i^{+}\equiv {|1⟩}_{\mathit{ii}}⟨0|$ and ${\sigma }_i^{-}\equiv {|0⟩}_{\mathit{ii}}⟨1|$. As before in the introduction, we denote the cavity photon annihilation operator by c. Moreover, $k_{\mathrm{cav}}$ is the wavenumber of the photons inside the resonator and ${\omega }_0$, ${\omega }_{\mathrm{cav}}$ and g denote the atomic frequency, the cavity frequency and the (real) atom-cavity coupling constant, respectively. In the interaction picture with respect to the Hamiltonian $H_0$ of the cavity field and the atomic states,(11) $$\begin{array}{cc}\hfill H_0& =ħ{\omega }_0{c}^{†}c+\sum _{i=1}^Nħ{\omega }_0{\sigma }_i^{+}{\sigma }_i^{-},\hfill \end{array}$$(11)

the total Hamiltonian of the system can now be written as(12) $$\begin{array}{cc}\hfill H_{\mathrm{I}}& =\sum _{i=1}^N\frac{ħ\mathrm{\Omega }}{2}{\sigma }_i^{-}+ħgsin(k_{\mathrm{cav}}{r}_i)c{\sigma }_i^{+}+\text{H.c.}\hfill \\ \hfill & +ħ\delta c^{†}c+H_{\mathrm{vib}}\hfill \end{array}$$(12)

in the usual dipole and rotating wave approximations and with the atom-cavity detuning $\delta$ defined as(13) $$\begin{array}{cc}\hfill \delta & ={\omega }_{\mathrm{cav}}-{\omega }_0.\hfill \end{array}$$(13)

To maximise the cooling process, this detuning should be equal to the phonon frequency $\delta =\nu$, as indicated in Equation (7(7) $$\begin{array}{cc}\hfill {\omega }_{\mathrm{cav}}-{\omega }_0& =\nu ,\hfill \end{array}$$(7) 13(13) $$\begin{array}{cc}\hfill \delta & ={\omega }_{\mathrm{cav}}-{\omega }_0.\hfill \end{array}$$(13) ) and as illustrated in Figure 3. For simplicity we neglect in the following direct dipole-dipole interactions between atomic particles. Since the applied laser field is in resonance with the atoms, transitions which do not involve the cavity field are not expected to change mean phonon numbers. Later, in Section 3.3.1, we approximate the electronic states of the atoms by a stationary state. Dipole-dipole interactions might slightly change this stationary state but would not change the predictions in Section 3.3.2 qualitatively.

Moreover, to maximise the coupling of the atomic motion to the cavity field mode, we assume in the following that the centre of the asymmetric trap R coincides with a node of the resonator field, as illustrated in Figure 2. This condition implies $sin(k_{\mathrm{cav}}R)=0$. As mentioned already above, we assume in addition that the trap localises the atoms within an area much smaller than the relevant wavelength of their electronic transition. Taking this and the fact that $r_i=R+x_i$ into account and applying the so-called Lamb-Dicke approximation yields(14) $$\begin{array}{cc}\hfill sin(k_{\mathrm{cav}}{r}_i)& =\eta (b_i+b_i^{†})\hfill \end{array}$$(14)

up to terms which scale as ${\eta }^3$, where $\eta$ is the so-called Lamb-Dicke parameter and much smaller than one. Hence(15) $$\begin{array}{cc}\hfill H_{\mathrm{I}}& =\sum _{i=1}^N\frac{ħ\mathrm{\Omega }}{2}{\sigma }_i^{-}+ħ\eta g(b_i+b_i^{†})c{\sigma }_i^{+}+\text{H.c.}\hfill \\ \hfill & +ħ\delta c^{†}c+H_{\mathrm{vib}}.\hfill \end{array}$$(15)

Finally, spontaneous emission from excited atoms and the leakage of photons through the cavity mirrors, is in the following taken into account by the quantum optical master equation [22,28](16) $$\begin{array}{cc}\hfill {\dot{\rho }}_{\mathrm{I}}& =-\frac{\mathrm{i}}{ħ}\left[H_{\mathrm{I}},{\rho }_{\mathrm{I}}\right]+\frac{\kappa }{2}\left(2c{\rho }_{\mathrm{I}}{c}^{†}-c^{†}c{\rho }_{\mathrm{I}}-{\rho }_{\mathrm{I}}{c}^{†}c\right)\hfill \\ \hfill & +\sum _{i=1}^N\frac{\mathrm{\Gamma }}{2}\left(2{\sigma }_i^{-}{\rho }_{\mathrm{I}}{\sigma }_i^{+}-{\sigma }_i^{+}{\sigma }_i^{-}{\rho }_{\mathrm{I}}-{\rho }_{\mathrm{I}}{\sigma }_i^{+}{\sigma }_i^{-}\right),\hfill \end{array}$$(16)

where $\mathrm{\Gamma }$ and $\kappa$ are the usual atom and cavity decay rates and where ${\rho }_{\mathrm{I}}$ denotes the density matrix of the atom-phonon-photon system in the interaction picture with respect to $H_0$ in Equation (11(11) $$\begin{array}{cc}\hfill H_0& =ħ{\omega }_0{c}^{†}c+\sum _{i=1}^Nħ{\omega }_0{\sigma }_i^{+}{\sigma }_i^{-},\hfill \end{array}$$(11) ).

3. Analysis of the cooling process

We now have all the theoretical tools needed to analyse the dynamics of the trapped atomic gas during the proposed cyclic two-stage cooling process. However before doing so, let us have a closer look at the most relevant expectation values of the experimental setup in Figure 2. Usually, in laser cooling, there is only one slowly-evolving expectation value, namely the mean phonon number m per particle (cf. Equation (3(3) $$\begin{array}{ccc}\hfill m& \equiv \hfill & \hfill \frac{1}{N}\sum _{i=1}^N⟨b_i^{†}{b}_i⟩,\\ \hfill \zeta & \equiv \hfill & \hfill \frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}⟨b_i^{†}{b}_j⟩.\end{array}$$(3) )) [2]. However, as we shall see below, the cooling scheme that we propose here has two slowly-evolving phonon expectation values, which play a crucial role in the cooling process. One is the mean phonon number m. The second variable is the many-particle phonon expectation value $\zeta$ which we introduced in Equation (3(3) $$\begin{array}{ccc}\hfill m& \equiv \hfill & \hfill \frac{1}{N}\sum _{i=1}^N⟨b_i^{†}{b}_i⟩,\\ \hfill \zeta & \equiv \hfill & \hfill \frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}⟨b_i^{†}{b}_j⟩.\end{array}$$(3) ). In this section, we analyse the dynamics of both variables using standard quantum optical approximations.

3.1. The physical meaning of m and $\zeta$

But before doing so, let us have a closer look at the physical meaning of m and $\zeta$ in Equation (3(3) $$\begin{array}{ccc}\hfill m& \equiv \hfill & \hfill \frac{1}{N}\sum _{i=1}^N⟨b_i^{†}{b}_i⟩,\\ \hfill \zeta & \equiv \hfill & \hfill \frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}⟨b_i^{†}{b}_j⟩.\end{array}$$(3) ). In this paper, we refer to m and $\zeta$ as phonon expectation values, since they are expectation values of the bosonic operators $b_i$ in Equation (9(9) $$\begin{array}{ccc}\hfill b_i& \equiv \hfill & \hfill \sqrt{\frac{M\nu }{2ħ}}{x}_i+\mathrm{i}\sqrt{\frac{1}{2Mħ\nu }}{p}_i,\\ \hfill b_i^{†}& \equiv \hfill & \hfill \sqrt{\frac{M\nu }{2ħ}}{x}_i-\mathrm{i}\sqrt{\frac{1}{2Mħ\nu }}{p}_i,\end{array}$$(9) ). If the particles were confined by a harmonic trapping potential, m would be a direct measure of their mean vibrational energy. However, in an asymmetric trap, the expectation value for the vibrational energy of the trapped atoms is given by $⟨H_{\mathrm{vib}}⟩$ and no longer coincides automatically with $Nħ\nu m$. Comparing Equations (3(3) $$\begin{array}{ccc}\hfill m& \equiv \hfill & \hfill \frac{1}{N}\sum _{i=1}^N⟨b_i^{†}{b}_i⟩,\\ \hfill \zeta & \equiv \hfill & \hfill \frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}⟨b_i^{†}{b}_j⟩.\end{array}$$(3) ) and (10(10) $$\begin{array}{cc}\hfill H_{\mathrm{vib}}& =\sum _{i=1}^Nħ\nu b_i^{†}{b}_i+ħ\mu {(b_i+b_i^{†})}^3,\hfill \end{array}$$(10) ), one can show that(17) $$\begin{array}{cc}\hfill ⟨H_{\mathrm{vib}}⟩& =Nħ\nu m+\sum _{i=1}^Nħ\mu 〈{(b_i+b_i^{†})}^3〉.\hfill \end{array}$$(17)

This shows that $⟨H_{\mathrm{vib}}⟩=Nħ\nu m$ only applies to a very good approximation when(18) $$\begin{array}{cc}\hfill \mu \ll \nu \&\text{and/or}& m\ll 1.\hfill \end{array}$$(18)

In the remainder of this paper, we restrict ourselves to this parameter regime and consider only slightly asymmetric trapping potentials with $\mu \ll \nu$. Moreover, at the end of the cooling process which we analyse in this paper predict final phonon numbers $m_{\mathrm{final}}$ which are much smaller than one. Hence we can safely assume here that the atoms are eventually cooled down to very low temperatures as long as m becomes eventually very small. In general, cooling only occurs as long as $⟨H_{\mathrm{vib}}⟩$ decreases in time.

Next let us have a closer look at the physical meaning of the phonon expectation value $\zeta$. Combining Equations (3(3) $$\begin{array}{ccc}\hfill m& \equiv \hfill & \hfill \frac{1}{N}\sum _{i=1}^N⟨b_i^{†}{b}_i⟩,\\ \hfill \zeta & \equiv \hfill & \hfill \frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}⟨b_i^{†}{b}_j⟩.\end{array}$$(3) ) and (9(9) $$\begin{array}{ccc}\hfill b_i& \equiv \hfill & \hfill \sqrt{\frac{M\nu }{2ħ}}{x}_i+\mathrm{i}\sqrt{\frac{1}{2Mħ\nu }}{p}_i,\\ \hfill b_i^{†}& \equiv \hfill & \hfill \sqrt{\frac{M\nu }{2ħ}}{x}_i-\mathrm{i}\sqrt{\frac{1}{2Mħ\nu }}{p}_i,\end{array}$$(9) ), one can show that(19) $$\begin{array}{cc}\hfill \zeta & =\frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}\frac{M\nu }{2ħ}⟨x_i{x}_j⟩+\frac{1}{2ħM\nu }⟨p_i{p}_j⟩\hfill \end{array}$$(19)

without any approximations. Assuming that all particles move independently and that the atomic confinement is so strong that there is no flux of particles out of the trap and applying semiclassical standard approximations, this implies that(20) $$\begin{array}{cc}\hfill ⟨x_i{x}_j⟩& =⟨x_i⟩⟨x_j⟩,\hfill \\ \hfill ⟨p_i{p}_j⟩& =⟨p_i⟩⟨p_j⟩=0.\hfill \end{array}$$(20)

Combining these two observations with Equation (19(19) $$\begin{array}{cc}\hfill \zeta & =\frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}\frac{M\nu }{2ħ}⟨x_i{x}_j⟩+\frac{1}{2ħM\nu }⟨p_i{p}_j⟩\hfill \end{array}$$(19) ) shows(21) $$\begin{array}{cc}\hfill \zeta & =\frac{M\nu }{2ħ}{⟨x⟩}^2\hfill \end{array}$$(21)

to a very good approximation, where $⟨x⟩$ denotes the mean distance of the trapped particles from the centre of the trap. When this equation applies, $\zeta$ is always positive and a measure for the distance of the atomic gas from the centre of the trap. For example, $⟨x⟩=0$ implies $\zeta =0$ and the accumulation of the atomic particles around the centre of the trap.

Having a closer look at Equation (19(19) $$\begin{array}{cc}\hfill \zeta & =\frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}\frac{M\nu }{2ħ}⟨x_i{x}_j⟩+\frac{1}{2ħM\nu }⟨p_i{p}_j⟩\hfill \end{array}$$(19) ), which is exact, we see that the classical analog of the phonon expectation value $\zeta$ is in general different from zero. For example, if the symmetry of the trap is broken and many more particles accumulate on one side of the centre of the trap than on the other, then the first term in Equation (19(19) $$\begin{array}{cc}\hfill \zeta & =\frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}\frac{M\nu }{2ħ}⟨x_i{x}_j⟩+\frac{1}{2ħM\nu }⟨p_i{p}_j⟩\hfill \end{array}$$(19) ) contains significantly more positive than negative contributions. If the particles move moreover in a synchronised fashion, such that the momenta of a majority of particles have the same sign, then the second term in Equation (19(19) $$\begin{array}{cc}\hfill \zeta & =\frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}\frac{M\nu }{2ħ}⟨x_i{x}_j⟩+\frac{1}{2ħM\nu }⟨p_i{p}_j⟩\hfill \end{array}$$(19) ) increases $\zeta$ even further. In the following, we ignore the possibility of synchronised motion and estimate $\zeta$ at the end of every displacement stage using Equation (21(21) $$\begin{array}{cc}\hfill \zeta & =\frac{M\nu }{2ħ}{⟨x⟩}^2\hfill \end{array}$$(21) ), thereby potentially underestimating the efficiency of the proposed cooling scheme. However, as we shall see below, at the end of every cooling stage, we have $\zeta =0$ which suggests the absence of synchronisation due to the common coupling of the particles to the cavity field. As one would expect for independently and freely moving particles, we then obtain a positive value for $\zeta$ at the end of each displacement stage which increases as the trap asymmetry increases.

Notice also that $\zeta$ can only be defined when there are at least two particles in the trap. For only one particle, Equation (21(21) $$\begin{array}{cc}\hfill \zeta & =\frac{M\nu }{2ħ}{⟨x⟩}^2\hfill \end{array}$$(21) ) does not apply and we always have $\zeta \equiv 0$ (cf. Equation (3(3) $$\begin{array}{ccc}\hfill m& \equiv \hfill & \hfill \frac{1}{N}\sum _{i=1}^N⟨b_i^{†}{b}_i⟩,\\ \hfill \zeta & \equiv \hfill & \hfill \frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}⟨b_i^{†}{b}_j⟩.\end{array}$$(3) )). This is why the cooling process, which we describe in this paper, is qualitatively different from single-particle laser-cooling schemes. Single-particle cavity cooling schemes cannot take advantage of a non-zero $\zeta$ as an additional resource.

3.2. Displacement stage

As we shall see below, at the end of every cooling stage, and therefore also at the beginning of every displacement stage, we have $\zeta =0$, while the average mean phonon number m is in general non-zero. The main purpose of the displacement stage is to increase $\zeta$, while keeping m essentially the same. In this subsection, we show that this task can be achieved by allowing the particles oscillate freely and independently within their asymmetric trapping potential. To do so, we use the conservation of the vibrational energy and Equation (21(21) $$\begin{array}{cc}\hfill \zeta & =\frac{M\nu }{2ħ}{⟨x⟩}^2\hfill \end{array}$$(21) ) to estimate the dynamics of $\zeta$ and m. In the absence of any laser driving and other quantum effects, the atomic particles move according to the laws of classical mechanics and there is nothing quantum about their motion.

Unfortunately, analysing the dynamics of particles inside an asymmetric trapping potential, even classically, is not straightforward. In this subsection, we therefore employ an approximation. To determine the average position x(t) of a single trapped particle after a sufficiently long time t, we assume in the following that(22) $$\begin{array}{cc}\hfill x(t)& =\frac{x_{\min }+x_{\max }}{2}\hfill \end{array}$$(22)

to a very good approximation. Here $x_{\min }$ and $x_{\max }$ are the return points of the free oscillations of the atom. These two positions can be calculated relatively easily, since they correspond to points of vanishing momentum p. Combining Equations (9(9) $$\begin{array}{ccc}\hfill b_i& \equiv \hfill & \hfill \sqrt{\frac{M\nu }{2ħ}}{x}_i+\mathrm{i}\sqrt{\frac{1}{2Mħ\nu }}{p}_i,\\ \hfill b_i^{†}& \equiv \hfill & \hfill \sqrt{\frac{M\nu }{2ħ}}{x}_i-\mathrm{i}\sqrt{\frac{1}{2Mħ\nu }}{p}_i,\end{array}$$(9) ) and (10(10) $$\begin{array}{cc}\hfill H_{\mathrm{vib}}& =\sum _{i=1}^Nħ\nu b_i^{†}{b}_i+ħ\mu {(b_i+b_i^{†})}^3,\hfill \end{array}$$(10) ) and setting $p=0$, we find that(23) $$\begin{array}{cc}\hfill H_{\mathrm{vib}}& =\frac{M{\nu }^2}{2}{x}^2+ħ\mu {\left(\frac{2M\nu }{ħ}\right)}^{3/2}{x}^3\hfill \end{array}$$(23)

for $x=x_{\min }$ and $x=x_{\max }$. Moreover, we know that the vibrational energy of the trapped particle $H_{\mathrm{vib}}$, which is conserved, equals(24) $$\begin{array}{cc}\hfill H_{\mathrm{vib}}(0)& =ħ\nu m(0)\hfill \end{array}$$(24)

at the beginning of the displacement stage, since $\zeta$ and therefore also x are both zero at $t=0$. Combining Equations (22(22) $$\begin{array}{cc}\hfill x(t)& =\frac{x_{\min }+x_{\max }}{2}\hfill \end{array}$$(22) )–(24(24) $$\begin{array}{cc}\hfill H_{\mathrm{vib}}(0)& =ħ\nu m(0)\hfill \end{array}$$(24) ) yields(25) $$\begin{array}{cc}\hfill x(t)& =-\frac{4\mu }{\nu }{\left(\frac{2ħ}{M\nu }\right)}^{1/2}m(0)\hfill \end{array}$$(25)

up to first order in $\mu /\nu$. Substituting this result into Equation (21(21) $$\begin{array}{cc}\hfill \zeta & =\frac{M\nu }{2ħ}{⟨x⟩}^2\hfill \end{array}$$(21) ), we finally see that(26) $$\begin{array}{cc}\hfill \zeta (t)& ={\left(\frac{4\mu }{\nu }\right)}^{2}m{(0)}^2\hfill \end{array}$$(26)

to a good approximation. As one would expect, $\zeta (t)$ is positive and non-zero at the end of every displacement stage, unless $H_{\mathrm{vib}}$ tends to zero. Moreover, $\zeta$ in Equation (26(26) $$\begin{array}{cc}\hfill \zeta (t)& ={\left(\frac{4\mu }{\nu }\right)}^{2}m{(0)}^2\hfill \end{array}$$(26) ) increases rapidly as the vibrational energy, i.e. m(0), of the atomic particles increases. The hotter the atoms are, the further away they accumulate on average from the centre of the trap.

Combining Equations (9(9) $$\begin{array}{ccc}\hfill b_i& \equiv \hfill & \hfill \sqrt{\frac{M\nu }{2ħ}}{x}_i+\mathrm{i}\sqrt{\frac{1}{2Mħ\nu }}{p}_i,\\ \hfill b_i^{†}& \equiv \hfill & \hfill \sqrt{\frac{M\nu }{2ħ}}{x}_i-\mathrm{i}\sqrt{\frac{1}{2Mħ\nu }}{p}_i,\end{array}$$(9) ), (10(10) $$\begin{array}{cc}\hfill H_{\mathrm{vib}}& =\sum _{i=1}^Nħ\nu b_i^{†}{b}_i+ħ\mu {(b_i+b_i^{†})}^3,\hfill \end{array}$$(10) ), (21(21) $$\begin{array}{cc}\hfill \zeta & =\frac{M\nu }{2ħ}{⟨x⟩}^2\hfill \end{array}$$(21) ) and (26(26) $$\begin{array}{cc}\hfill \zeta (t)& ={\left(\frac{4\mu }{\nu }\right)}^{2}m{(0)}^2\hfill \end{array}$$(26) ), we find that the vibrational energy of a single atomic particle equals(27) $$\begin{array}{cc}\hfill H_{\mathrm{vib}}(t)& =ħ\nu m(t)+8ħ\mu \zeta {(t)}^{3/2}\hfill \end{array}$$(27)

at the end of the displacement stage. Conservation of the vibrational energy during the free evolution of the trapped particle implies that this expression is the same as $H_{\mathrm{vib}}(0)$ in Equation (24(24) $$\begin{array}{cc}\hfill H_{\mathrm{vib}}(0)& =ħ\nu m(0)\hfill \end{array}$$(24) ). Using this and Equation (26(26) $$\begin{array}{cc}\hfill \zeta (t)& ={\left(\frac{4\mu }{\nu }\right)}^{2}m{(0)}^2\hfill \end{array}$$(26) ), one can now show that(28) $$\begin{array}{cc}\hfill m(t)& =m(0)\left[1-2{\left(\frac{4\mu }{\nu }\right)}^{4}m{(0)}^2\right]\hfill \end{array}$$(28)

to a very good approximation. For sufficiently small trap asymmetries $\mu /\nu$, the average mean phonon number m remains essentially the same.

3.3. Cooling stage

Next we analyse the dynamics of the expectation values m and $\zeta$ when the cooling laser is turned on. To do so, we first simplify the full master equation in Equation (16(16) $$\begin{array}{cc}\hfill {\dot{\rho }}_{\mathrm{I}}& =-\frac{\mathrm{i}}{ħ}\left[H_{\mathrm{I}},{\rho }_{\mathrm{I}}\right]+\frac{\kappa }{2}\left(2c{\rho }_{\mathrm{I}}{c}^{†}-c^{†}c{\rho }_{\mathrm{I}}-{\rho }_{\mathrm{I}}{c}^{†}c\right)\hfill \\ \hfill & +\sum _{i=1}^N\frac{\mathrm{\Gamma }}{2}\left(2{\sigma }_i^{-}{\rho }_{\mathrm{I}}{\sigma }_i^{+}-{\sigma }_i^{+}{\sigma }_i^{-}{\rho }_{\mathrm{I}}-{\rho }_{\mathrm{I}}{\sigma }_i^{+}{\sigma }_i^{-}\right),\hfill \end{array}$$(16) ) to an effective master equation that has many similarities with the master equation of the idealised toy model in Section 1. Afterwards we analyse its effective time evolution. It is shown that the mean number of phonons per particle m drops very rapidly by an amount which equals the initial value of the phonon expectation value $\zeta$. In other words, the dynamics of the experimental setup are analogous to the dynamics described in Section 1, even when off-resonant interaction terms in the system Hamiltonian, additional decay channels and the asymmetry of the trap are taken into account. What we find is that the change of the mean phonon number m depends only on the sign of $\zeta$. As long as $\zeta$ is positive, we predict a reduction of m with a collectively enhanced cooling rate, as described in Equation (6(6) $$\begin{array}{cc}\hfill m(t)& =m(0)-\zeta (0),\hfill \\ \hfill \zeta (t)& =0\hfill \end{array}$$(6) ).

3.3.1. An effective master equation

As illustrated in Figure 3, the role of the atoms is to bridge the energy gap between phonons and cavity photons, thereby helping to create an effective resonant direct exchange interaction between these different types of particles. Given the parameters of typical atom-cavity experiments, where¹(29) $$\begin{array}{cc}\hfill \mathrm{\Omega },\mathrm{\Gamma }& \gg \nu ,\kappa ,\delta ,g,\hfill \end{array}$$(29)

the conversion is usually relatively inefficient, even when it is collectively enhanced, while the electronic states of the atoms evolve on a relatively fast time scale. In the following, we take this into account and approximate the state of the atoms by the stationary state, which they would assume in the absence of any off-resonant coupling to the cavity field and phonon modes. More concretely, we assume in the following that the atomic density matrix ${\rho }_{\mathrm{A}}(t)$ evolves in the interaction picture according to the master equation(30) $$\begin{array}{cc}\hfill {\dot{\rho }}_{\mathrm{A}}& =-\frac{\mathrm{i}\mathrm{\Omega }}{2}\left[{\sigma }_i^{-}+{\sigma }_i^{-},{\rho }_{\mathrm{A}}\right]\hfill \\ \hfill & +\sum _{i=1}^N\frac{\mathrm{\Gamma }}{2}\left(2{\sigma }_i^{-}{\rho }_{\mathrm{A}}{\sigma }_i^{+}-{\sigma }_i^{+}{\sigma }_i^{-}{\rho }_{\mathrm{A}}-{\rho }_{\mathrm{A}}{\sigma }_i^{+}{\sigma }_i^{-}\right)\hfill \end{array}$$(30)

and approximate the state of the atoms by the stationary state ${\rho }_{\mathrm{A}}^{\mathrm{ss}}$ with ${\dot{\rho }}_{\mathrm{A}}=0$.

In other words, we assume in the following that the density matrix ${\rho }_{\mathrm{I}}(t)$ can always be written as(31) $$\begin{array}{cc}\hfill {\rho }_{\mathrm{I}}(t)& ={\rho }_{\mathrm{BC}}(t)\otimes {\rho }_{\mathrm{A}}^{\mathrm{ss}},\hfill \end{array}$$(31)

where ${\rho }_{\mathrm{BC}}(t)$ denotes the density matrix of phonons and cavity photons at time t. If the atoms reach their stationary state ${\rho }_{\mathrm{A}}^{\mathrm{ss}}$ very rapidly without changing the statistical properties of the phonons and photons in the system, the master Equation (16(16) $$\begin{array}{cc}\hfill {\dot{\rho }}_{\mathrm{I}}& =-\frac{\mathrm{i}}{ħ}\left[H_{\mathrm{I}},{\rho }_{\mathrm{I}}\right]+\frac{\kappa }{2}\left(2c{\rho }_{\mathrm{I}}{c}^{†}-c^{†}c{\rho }_{\mathrm{I}}-{\rho }_{\mathrm{I}}{c}^{†}c\right)\hfill \\ \hfill & +\sum _{i=1}^N\frac{\mathrm{\Gamma }}{2}\left(2{\sigma }_i^{-}{\rho }_{\mathrm{I}}{\sigma }_i^{+}-{\sigma }_i^{+}{\sigma }_i^{-}{\rho }_{\mathrm{I}}-{\rho }_{\mathrm{I}}{\sigma }_i^{+}{\sigma }_i^{-}\right),\hfill \end{array}$$(16) ) can be used to show that, for a small time interval $\mathrm{\Delta }t$,(32) $$\begin{array}{cc}\hfill {\dot{\rho }}_{\mathrm{BC}}(t+\mathrm{\Delta }t)& ={\mathrm{Tr}}_{\mathrm{A}}\left({\rho }_{\mathrm{I}}(t)+{\dot{\rho }}_{\mathrm{I}}(t)\mathrm{\Delta }t\right),\hfill \end{array}$$(32)

where ${\mathrm{Tr}}_{\mathrm{A}}(·)$ denotes the trace operation over the state space of the electronic states of the atoms. Combining this equation with Equations (16(16) $$\begin{array}{cc}\hfill {\dot{\rho }}_{\mathrm{I}}& =-\frac{\mathrm{i}}{ħ}\left[H_{\mathrm{I}},{\rho }_{\mathrm{I}}\right]+\frac{\kappa }{2}\left(2c{\rho }_{\mathrm{I}}{c}^{†}-c^{†}c{\rho }_{\mathrm{I}}-{\rho }_{\mathrm{I}}{c}^{†}c\right)\hfill \\ \hfill & +\sum _{i=1}^N\frac{\mathrm{\Gamma }}{2}\left(2{\sigma }_i^{-}{\rho }_{\mathrm{I}}{\sigma }_i^{+}-{\sigma }_i^{+}{\sigma }_i^{-}{\rho }_{\mathrm{I}}-{\rho }_{\mathrm{I}}{\sigma }_i^{+}{\sigma }_i^{-}\right),\hfill \end{array}$$(16) ) and (30(30) $$\begin{array}{cc}\hfill {\dot{\rho }}_{\mathrm{A}}& =-\frac{\mathrm{i}\mathrm{\Omega }}{2}\left[{\sigma }_i^{-}+{\sigma }_i^{-},{\rho }_{\mathrm{A}}\right]\hfill \\ \hfill & +\sum _{i=1}^N\frac{\mathrm{\Gamma }}{2}\left(2{\sigma }_i^{-}{\rho }_{\mathrm{A}}{\sigma }_i^{+}-{\sigma }_i^{+}{\sigma }_i^{-}{\rho }_{\mathrm{A}}-{\rho }_{\mathrm{A}}{\sigma }_i^{+}{\sigma }_i^{-}\right)\hfill \end{array}$$(30) ), we find that the time derivative ${\rho }_{\mathrm{BC}}(t)$ equals(33) $$\begin{array}{cc}\hfill {\dot{\rho }}_{\mathrm{BC}}(t)& ={\mathrm{Tr}}_{\mathrm{A}}\left({\dot{\rho }}_{\mathrm{I}}(t)\right).\hfill \end{array}$$(33)

Hence the density matrix ${\rho }_{\mathrm{BC}}$ evolves effectively according to the master equation(34) $$\begin{array}{cc}\hfill {\dot{\rho }}_{\mathrm{BC}}& ={\mathrm{Tr}}_{\mathrm{A}}\left(-\frac{\mathrm{i}}{ħ}\left[H_{\mathrm{I}},{\rho }_{\mathrm{BC}}\otimes {\rho }_{\mathrm{A}}^{\mathrm{ss}}\right]\right)\hfill \\ \hfill & +\frac{\kappa }{2}\left(2c{\rho }_{\mathrm{BC}}{c}^{†}-c^{†}c{\rho }_{\mathrm{BC}}-{\rho }_{\mathrm{BC}}{c}^{†}c\right).\hfill \end{array}$$(34)

This equation is of exactly the same form as Equation (2(2) $$\begin{array}{c}\hfill {\dot{\rho }}_{\mathrm{I}}=-\frac{\mathrm{i}}{ħ}\left[H_{\mathrm{ideal}},{\rho }_{\mathrm{I}}\right]+\frac{\kappa }{2}\left(2c{\rho }_{\mathrm{I}}{c}^{†}-c^{†}c{\rho }_{\mathrm{I}}-{\rho }_{\mathrm{I}}{c}^{†}c\right).\end{array}$$(2) ) but, in contrast to Equation (1(1) $$\begin{array}{cc}\hfill H_{\mathrm{ideal}}& =\sum _{i=1}^NħJ\left[b_i{c}^{†}+b_i^{†}c\right]\hfill \end{array}$$(1) ), the effective phonon-photon interaction Hamiltonian $H_{\mathrm{eff}}$ is now given by(35) $$\begin{array}{cc}\hfill H_{\mathrm{eff}}& =\sum _{i=1}^Nħ\eta g{}_i⟨0|{\rho }_{\mathrm{A}}^{\mathrm{ss}}{|1⟩}_i\left(b_i+b_i^{†}\right)c+\text{H.c.}\hfill \\ \hfill & +ħ\delta c^{†}c+H_{\mathrm{vib}},\hfill \end{array}$$(35)

where ${|0⟩}_i$ denotes the ground state of atom i, while ${|1⟩}_i$ denotes its excited state. Proceeding as described in Appendix B, analysing the atomic dynamics and calculating the corresponding stationary state in the presence of resonant laser driving and spontaneous photon emission, we finally find that(36) $$\begin{array}{cc}\hfill H_{\mathrm{eff}}& =\sum _{i=1}^N\mathrm{i}ħ\eta g_{\mathrm{eff}}\left(b_i+b_i^{†}\right)\left(c-c^{†}\right)+ħ\nu b_i^{†}{b}_i\hfill \\ \hfill & +ħ\mu {\left(b_i+b_i^{†}\right)}^3+ħ\delta c^{†}c\hfill \end{array}$$(36)

with the effective coupling constant $g_{\mathrm{eff}}$ given by(37) $$\begin{array}{cc}\hfill g_{\mathrm{eff}}& =-\frac{g\mathrm{\Gamma }\mathrm{\Omega }}{{\mathrm{\Gamma }}^2+2{\mathrm{\Omega }}^2}.\hfill \end{array}$$(37)

This Hamiltonian has many similarities with the toy model Hamiltonian $H_{\mathrm{eff}}$ in Equation (1(1) $$\begin{array}{cc}\hfill H_{\mathrm{ideal}}& =\sum _{i=1}^NħJ\left[b_i{c}^{†}+b_i^{†}c\right]\hfill \end{array}$$(1) ). Moreover, the master equation of the phonon and cavity photon system for the experimental setup in Figure 2 is the same as the toy model master equation in Equation (2(2) $$\begin{array}{c}\hfill {\dot{\rho }}_{\mathrm{I}}=-\frac{\mathrm{i}}{ħ}\left[H_{\mathrm{ideal}},{\rho }_{\mathrm{I}}\right]+\frac{\kappa }{2}\left(2c{\rho }_{\mathrm{I}}{c}^{†}-c^{†}c{\rho }_{\mathrm{I}}-{\rho }_{\mathrm{I}}{c}^{†}c\right).\end{array}$$(2) ). As we shall see below, for positive atom-cavity detunings $\delta$, the phonon dynamics of the proposed cooling scheme is essentially the same as the phonon dynamics described in Section 1. To show that this is indeed the case, we now analyse the dynamics of m and $\zeta$ with all additional coupling terms, detunings and decay channels taken into account.

3.3.2. Effective cooling equations

First we have a closer look at the dynamics of the slowly evolving variables m an $\zeta$. Using the above derived effective system dynamics (cf. Equations (32(32) $$\begin{array}{cc}\hfill {\dot{\rho }}_{\mathrm{BC}}(t+\mathrm{\Delta }t)& ={\mathrm{Tr}}_{\mathrm{A}}\left({\rho }_{\mathrm{I}}(t)+{\dot{\rho }}_{\mathrm{I}}(t)\mathrm{\Delta }t\right),\hfill \end{array}$$(32) 34(34) $$\begin{array}{cc}\hfill {\dot{\rho }}_{\mathrm{BC}}& ={\mathrm{Tr}}_{\mathrm{A}}\left(-\frac{\mathrm{i}}{ħ}\left[H_{\mathrm{I}},{\rho }_{\mathrm{BC}}\otimes {\rho }_{\mathrm{A}}^{\mathrm{ss}}\right]\right)\hfill \\ \hfill & +\frac{\kappa }{2}\left(2c{\rho }_{\mathrm{BC}}{c}^{†}-c^{†}c{\rho }_{\mathrm{BC}}-{\rho }_{\mathrm{BC}}{c}^{†}c\right).\hfill \end{array}$$(34) ) and (36(36) $$\begin{array}{cc}\hfill H_{\mathrm{eff}}& =\sum _{i=1}^N\mathrm{i}ħ\eta g_{\mathrm{eff}}\left(b_i+b_i^{†}\right)\left(c-c^{†}\right)+ħ\nu b_i^{†}{b}_i\hfill \\ \hfill & +ħ\mu {\left(b_i+b_i^{†}\right)}^3+ħ\delta c^{†}c\hfill \end{array}$$(36) )) and proceeding as described in Appendix B, one can now show that(38) $$\begin{array}{cc}\hfill \dot{m}& =\eta g_{\mathrm{eff}}{x}_{33}+3\mu y_{232},\hfill \\ \hfill \dot{\zeta }& =\eta g_{\mathrm{eff}}{x}_{33}+3\mu {\tilde{y}}_{223}\hfill \end{array}$$(38)

without any approximations. The definition of the x, y and $\tilde{y}$ variables on the right hand side of this equation can be found in Equation (B9(B9) $$\begin{array}{ccc}\hfill x_{\mathit{ab}}& \equiv \hfill & \hfill \frac{1}{N}\sum _{i=1}^N〈B_a^{(i)}{C}_b〉,\\ \hfill y_{\mathit{abc}}& \equiv \hfill & \hfill \frac{1}{N}\sum _{i=1}^N〈B_a^{(i)}{B}_b^{(i)}{B}_c^{(i)}〉,\\ \hfill {\tilde{y}}_{\mathit{abc}}& \equiv \hfill & \hfill \frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}〈B_a^{(i)}{B}_b^{(i)}{B}_c^{(j)}〉,\\ \hfill {\tilde{z}}_{\mathit{ab}}& \equiv \hfill & \hfill \frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}〈B_a^{(i)}{B}_b^{(j)}〉.\end{array}$$(B9) ) in Appendix B. Here $x_{33}$ and $y_{232}$ are single-particle expectation values, while ${\tilde{y}}_{223}$ denotes a many-particle phonon expectation value. As shown in Appendix B, the variables $x_{33}$, $y_{232}$ an ${\tilde{y}}_{223}$ evolve on a much faster time scale than m and $\zeta$ when(39) $$\begin{array}{cc}\hfill N\eta g_{\mathrm{eff}},\nu ,\kappa ,\delta & \gg \eta g_{\mathrm{eff}},\mu .\hfill \end{array}$$(39)

Taking this into account, proceeding as described in Appendix B and eliminating these variables adiabatically from the system dynamics, yields a closed set of two effective cooling equations. In analogy to Equation (4(4) $$\begin{array}{cc}\hfill \dot{m}=\dot{\zeta }& =-\frac{4N{J}^2}{\kappa }\zeta \hfill \end{array}$$(4) ), we find that(40) $$\begin{array}{cc}\hfill \dot{m}=\dot{\zeta }& =-A_N\zeta \hfill \end{array}$$(40)

with the collective rate $A_N$ given by(41) $$\begin{array}{c}\hfill A_N=\frac{N{(8\eta g\mathrm{\Gamma }\mathrm{\Omega })}^2\kappa \delta \nu }{\left[{\kappa }^4+16{({\delta }^2-{\nu }^2)}^2+8{\kappa }^2({\delta }^2+{\nu }^2)\right]{({\mathrm{\Gamma }}^2+2{\mathrm{\Omega }}^2)}^2}.\end{array}$$(41)

The factor N in Equation (41(41) $$\begin{array}{c}\hfill A_N=\frac{N{(8\eta g\mathrm{\Gamma }\mathrm{\Omega })}^2\kappa \delta \nu }{\left[{\kappa }^4+16{({\delta }^2-{\nu }^2)}^2+8{\kappa }^2({\delta }^2+{\nu }^2)\right]{({\mathrm{\Gamma }}^2+2{\mathrm{\Omega }}^2)}^2}.\end{array}$$(41) ) shows that the time evolution of the mean phonon number m is collectively enhanced. The system succeeds in taking advantage of the simultaneous coupling the atoms to the cavity field.

The effective rate equations in Equation (40(40) $$\begin{array}{cc}\hfill \dot{m}=\dot{\zeta }& =-A_N\zeta \hfill \end{array}$$(40) ) are formally the same as the rate equations in Equation (4(4) $$\begin{array}{cc}\hfill \dot{m}=\dot{\zeta }& =-\frac{4N{J}^2}{\kappa }\zeta \hfill \end{array}$$(4) ) and can be solved easily analogously. Suppose the atom-cavity detuning $\delta$ and the initial phonon expectation value $\zeta (0)$ are both positive and the cooling laser is turned on at $t=0$ for a time $t\gg 1/A_N$. Then $A_N$ is positive and m(t) and $\zeta (t)$ are again to a very good approximation given by Equation (6(6) $$\begin{array}{cc}\hfill m(t)& =m(0)-\zeta (0),\hfill \\ \hfill \zeta (t)& =0\hfill \end{array}$$(6) ). The cooling laser reduces the mean phonon number m by $\zeta (0)$. Playing the role of a resource, the variable $\zeta$ itself vanishes in the process which implies $⟨x(t)⟩=0$ at the end of the cooling stage (cf. Equation (21(21) $$\begin{array}{cc}\hfill \zeta & =\frac{M\nu }{2ħ}{⟨x⟩}^2\hfill \end{array}$$(21) )). Under the action of the laser cooling, the atomic gas returns on average to the centre of the trap.

One can easily see from Equation (41(41) $$\begin{array}{c}\hfill A_N=\frac{N{(8\eta g\mathrm{\Gamma }\mathrm{\Omega })}^2\kappa \delta \nu }{\left[{\kappa }^4+16{({\delta }^2-{\nu }^2)}^2+8{\kappa }^2({\delta }^2+{\nu }^2)\right]{({\mathrm{\Gamma }}^2+2{\mathrm{\Omega }}^2)}^2}.\end{array}$$(41) ) that the cooling laser reduces the mean phonon number m most rapidly when $\delta =\nu$. In this case, the effective cooling rate $A_N$ simplifies to(42) $$\begin{array}{cc}\hfill A_N& =\frac{N{(8\eta g\mathrm{\Gamma }\mathrm{\Omega }\nu )}^2}{\kappa ({\kappa }^2+16{\nu }^2){({\mathrm{\Gamma }}^2+2{\mathrm{\Omega }}^2)}^2}.\hfill \end{array}$$(42)

This is not unexpected. When $\delta =\nu$, the most likely transition in the dynamics of the system is the de-excitation of a particle accompanied by the annihilation of a phonon and the creation of a cavity photon [28]. Once the photon leaks out of the resonator, the phonon is permanently lost. Moreover, heating transitions, which involve the simultaneous excitation of the cavity and the creation of a phonon, are most unlikely in this case. However these heating and cooling transitions only affect the mean number of phonons m on a time scale much longer than $1/A_N$. For the cavity-mediated collective cooling scheme, which we propose here, they remain negligible. Instead, the dynamics of the cooling process is dominated by the collective effects which we illustrated in the Introduction with the help of a relatively simple toy model.

3.3.3. Some physical intuition

In Section 2.2, we assume that the atoms are strongly confined and that the centre of their trapping potential coincides with a node of the cavity field (cf. Figure 2). In this case, the above described cooling process can be understood relatively easily also intuitively. As illustrated in Figure 3, the applied laser field drives the atoms resonantly. Changing m therefore requires off-resonant interactions with the cavity field. For $\zeta =0$, the trapped particles experience equal amounts of cooling and heating, since their atom-cavity interaction becomes zero, and m remains constant in this case. However, for $\zeta \ne 0$ the dynamics of the system during cooling stages becomes less trivial. Cooling transitions move the atoms closer to the centre of the trap and reduce $\zeta$. Analogously, heating transitions move the atoms further away and increases $\zeta$. The dynamics of hotter particles is hence more intense. Eventually, all particles reach the centre of the trap which corresponds to a stable fixed point of their dynamics and where $\zeta =0$. The overall effect of the transfer of the atoms, from their initial position to the centre of the trap, is a loss of vibrational energy (cf. Equation (6(6) $$\begin{array}{cc}\hfill m(t)& =m(0)-\zeta (0),\hfill \\ \hfill \zeta (t)& =0\hfill \end{array}$$(6) )). To re-initiate the cooling process, the atoms need to be moved again away from the cavity node. Another displacement stage is needed.

4. The final mean phonon number

To achieve very low temperatures, we propose to alternate cooling stages with displacement stages as illustrated in Figure 1. During each displacement stage, the non-interacting atomic gas evolves freely and naturally accumulates a non-zero phonon coherence $\zeta$ without increasing its mean phonon number m. During each cooling stage, $\zeta$ tends to zero, while m decreases by the same amount in the process (cf. Equation (6(6) $$\begin{array}{cc}\hfill m(t)& =m(0)-\zeta (0),\hfill \\ \hfill \zeta (t)& =0\hfill \end{array}$$(6) )). The hotter the atoms are, the larger their initial phonon coherence $\zeta$ (cf. Equation (26(26) $$\begin{array}{cc}\hfill \zeta (t)& ={\left(\frac{4\mu }{\nu }\right)}^{2}m{(0)}^2\hfill \end{array}$$(26) )). Eventually, both expectations values, m and $\zeta$, become too small to dominate the dynamics of the system during cooling stages (cf. Equation (40(40) $$\begin{array}{cc}\hfill \dot{m}=\dot{\zeta }& =-A_N\zeta \hfill \end{array}$$(40) )). At this point, other cooling and heating processes, which occur only on a time scale given by ${\eta }^2$, are no longer negligible. These transitions have previously been analysed in detail in Ref. [28] for example.

To get at least an estimate for the final mean phonon number $m_{\mathrm{final}}$ of the proposed cooling process, we now combine Equations (26(26) $$\begin{array}{cc}\hfill \zeta (t)& ={\left(\frac{4\mu }{\nu }\right)}^{2}m{(0)}^2\hfill \end{array}$$(26) ) and (40(40) $$\begin{array}{cc}\hfill \dot{m}=\dot{\zeta }& =-A_N\zeta \hfill \end{array}$$(40) ) with Equation (39) in Ref. [28]. Doing so we find that the time derivative of m at the beginning of the next cooling stage equals(43) $$\begin{array}{cc}\hfill \dot{m}& =-\frac{16{\mu }^2{A}_N}{{\nu }^2}{m}^2-{\gamma }_{c}m+c\hfill \end{array}$$(43)

also to a very good approximation, where ${\gamma }_c$ and c are two constants of order ${\eta }^2$. Assuming that the cooling process ends, when $\dot{m}=0$, we find that(44) $$\begin{array}{cc}\hfill m_{\mathrm{final}}& =\frac{\nu }{4\mu }{\left(\frac{c}{A_N}\right)}^{1/2}\hfill \end{array}$$(44)

to a very good approximation. Most importantly, the right hand side of this equation scales as $1/N^{1/2}$ and tends to zero in the many-particle limit. For example, if we choose $\mu /\nu =0.01$ and $N={10}^6$, we have $\nu /(4\mu N^{1/2})=0.025$ which is much smaller than one. The cooling process, which we propose in this paper, therefore has the potential to cool an atomic gas to much lower temperatures than other laser cooling schemes.

For more traditional cavity-mediated laser cooling schemes, the final temperature of the atoms scales essentially as ${(\kappa /\nu )}^2$, which is independent of the number N of atoms inside the trap. For example, substituting c from Equation (41) in Ref. [28] into Equation (44(44) $$\begin{array}{cc}\hfill m_{\mathrm{final}}& =\frac{\nu }{4\mu }{\left(\frac{c}{A_N}\right)}^{1/2}\hfill \end{array}$$(44) ) yields(45) $$\begin{array}{cc}\hfill m_{\mathrm{final}}& =\frac{\kappa }{16\mu N^{1/2}}\hfill \end{array}$$(45)

when $\delta =\nu$ and $\mathrm{\Omega },\nu \ll \kappa ,\mathrm{\Gamma }$. As one would expect, collective cooling of an atomic gas to very low temperatures with the cavity-mediated laser cooling scheme which we propose here requires a not-too-small asymmetry of the trapping potential. Moreover the particles need to be strongly confined and the cavity decay rate $\kappa$ should be relatively small. Although the proposed cooling scheme might not require a cavity, which operates in the strong-coupling regime, the cavity should not be too bad either, as one would expect.

5. Conclusions

This paper proposes a cyclic two-stage laser cooling scheme to transfer a non-interacting atomic gas to very low temperatures. The corresponding experimental setup, in which the atoms are confined to an asymmetric trapping potential inside an optical resonator, is shown in Figure 2. The role of the cavity and the external laser field is to remove vibrational energy from the atoms and to release it in the form of photons into the environment. The role of the trap asymmetry is to create a positive phonon coherence $\zeta$ between laser interactions without affecting the free energy of the trapped particles. Our calculations assume that the atoms are strongly confined within a wavelength of the resonator. Moreover we assume that the asymmetry of the trap can be described by a cubic term with a constant $\mu$, which is much smaller than the phonon frequency $\nu$. Although the parameter regime, which we consider in this paper, is relatively standard for cavity-mediated laser-cooling (cf. Equations (18(18) $$\begin{array}{cc}\hfill \mu \ll \nu \&\text{and/or}& m\ll 1.\hfill \end{array}$$(18) ), (29(29) $$\begin{array}{cc}\hfill \mathrm{\Omega },\mathrm{\Gamma }& \gg \nu ,\kappa ,\delta ,g,\hfill \end{array}$$(29) ) and (39(39) $$\begin{array}{cc}\hfill N\eta g_{\mathrm{eff}},\nu ,\kappa ,\delta & \gg \eta g_{\mathrm{eff}},\mu .\hfill \end{array}$$(39) )), we find that the usual temperature limits do not apply. Instead, the final vibrational energy of each atom scales as one over the square root of the total number of particles N (cf. Equation (45(45) $$\begin{array}{cc}\hfill m_{\mathrm{final}}& =\frac{\kappa }{16\mu N^{1/2}}\hfill \end{array}$$(45) )) which vanishes in the infinitely-many particle limit.

One drawback of the proposed cooling scheme is that its cools the particles in only one dimension. However, as long as there is an energy exchange between different directions of motion, draining kinetic energy out of the x direction eventually also reduces the motion of the atoms in the y and in the z direction. Another drawback of the proposed cooling scheme is that experimental parameters need to be chosen relatively carefully. Equation (45(45) $$\begin{array}{cc}\hfill m_{\mathrm{final}}& =\frac{\kappa }{16\mu N^{1/2}}\hfill \end{array}$$(45) ) shows that the particle number N should be relatively large, while the cavity decay rate $\kappa$ and the trap asymmetry $\mu$ should not be too large and too small, respectively, in order to outperform standard laser cooling schemes. Nevertheless, the non-destructive cooling mechanism, which we propose in this paper, could become an important tool for quantum technological applications.

Acknowledgements

AB would like to dedicate this paper to her dear and very highly-regarded colleague, Danny Segal, who she met and interacted with during her time at Imperial College London. We also acknowledge inspiring discussions with P. Grangier and thank him for his constructive feedback on the manuscript. Statement of compliance with EPSRC policy framework on research data: This publication is theoretical work that does not require supporting research data.

Additional information

Funding

This work was supported by the UK Engineering and Physical Sciences Research Council EPSRC [grant ref. EP/H048901/1]; British Council and DST Government of India.

Notes

No potential conflict of interest was reported by the authors.

1 The collective cavity cooling scheme which we propose here also operates in the bad-cavity limit where g, $\kappa$ and $\mathrm{\Gamma }$ are all of similar size [33]. However making the assumption below simplifies the following calculations and helps us to point out its underlying mechanisms more clearly.

Appendix 1

Derivation of effective cooling equations for a ideal model

Using Equation (2(2) $$\begin{array}{c}\hfill {\dot{\rho }}_{\mathrm{I}}=-\frac{\mathrm{i}}{ħ}\left[H_{\mathrm{ideal}},{\rho }_{\mathrm{I}}\right]+\frac{\kappa }{2}\left(2c{\rho }_{\mathrm{I}}{c}^{†}-c^{†}c{\rho }_{\mathrm{I}}-{\rho }_{\mathrm{I}}{c}^{†}c\right).\end{array}$$(2) ), one can show that the time derivative of the expectation value $⟨A⟩$ of an observable A of our idealised cooling scheme equals(A1) $$\begin{array}{cc}\hfill ⟨\dot{A}⟩& =-\frac{\mathrm{i}}{ħ}〈\left[A,H_{\mathrm{ideal}}\right]〉+\frac{\kappa }{2}〈[c^{†},A]c+c^{†}[A,c]〉.\hfill \end{array}$$(A1)

For the parameter regime given in Equation (5(5) $$\begin{array}{c}\hfill \kappa \gg J\mathrm{and}N\gg 1.\end{array}$$(5) ), the dynamics of m and $\zeta$ is therefore to a very good approximation given by the closed set of rate equations(A2) $$\begin{array}{cc}\hfill \dot{m}=\dot{\zeta }& =Jx,\hfill \\ \hfill \dot{x}& =-2NJ\zeta -\frac{\kappa }{2}x\hfill \end{array}$$(A2)

with the additional multi-particle quantum coherence x defined as(A3) $$\begin{array}{ccc}\hfill x& \equiv \hfill & \hfill \frac{\mathrm{i}}{N}\sum _{i=1}^N⟨b_i{c}^{†}-b_i^{†}c⟩.\end{array}$$(A3)

The only approximation applied here is to neglect terms which scale as J in the time derivative of x. Since $N\gg 1$, the variable x evolves on a much faster time scale than the expectation values m and $\zeta$. Hence x can be adiabatically eliminated from the system dynamics. Doing so and assuming $\dot{x}=0$, we find that(A4) $$\begin{array}{cc}\hfill x& =-\frac{4NJ}{\kappa }\zeta \hfill \end{array}$$(A4)

to a very good approximation. Substituting this result into the differential equations for m and $\zeta$ in Equation (A2(A2) $$\begin{array}{cc}\hfill \dot{m}=\dot{\zeta }& =Jx,\hfill \\ \hfill \dot{x}& =-2NJ\zeta -\frac{\kappa }{2}x\hfill \end{array}$$(A2) ) eventually yields the closed set of two effective cooling equations in Equation (4(4) $$\begin{array}{cc}\hfill \dot{m}=\dot{\zeta }& =-\frac{4N{J}^2}{\kappa }\zeta \hfill \end{array}$$(4) ).

Appendix 2

Derivation of effective cooling equations for the proposed cooling scheme

As pointed out already in Section 3.3, the purpose of this appendix is to derive a closed set of rate equations, for the time evolution of the expectation values m and $\zeta$ which we introduced in Equation (19(19) $$\begin{array}{cc}\hfill \zeta & =\frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}\frac{M\nu }{2ħ}⟨x_i{x}_j⟩+\frac{1}{2ħM\nu }⟨p_i{p}_j⟩\hfill \end{array}$$(19) ). To do so, we first analyse the zeroth order dynamics of the atoms in the presence of resonant laser driving and spontaneous photon emission. Afterwards, we introduce further expectation values which we need to consider in addition to the slowly evolving expectation values $\zeta$ an m. Afterwards, we have a closer look at the time derivatives of these relatively rapidly evolving variables and calculate them with the help of an adiabatic elimination technique which is routinely used in quantum optics. The following results apply to a very goo approximation, as long as the system parameters obey the inequalities summarised in Equation (29(29) $$\begin{array}{cc}\hfill \mathrm{\Omega },\mathrm{\Gamma }& \gg \nu ,\kappa ,\delta ,g,\hfill \end{array}$$(29) ).

Atomic dynamics

First, we calculate the expectation value ${}_i⟨0|{\rho }_{\mathrm{A}}^{\mathrm{ss}}{|1⟩}_i$ in Equation (35(35) $$\begin{array}{cc}\hfill H_{\mathrm{eff}}& =\sum _{i=1}^Nħ\eta g{}_i⟨0|{\rho }_{\mathrm{A}}^{\mathrm{ss}}{|1⟩}_i\left(b_i+b_i^{†}\right)c+\text{H.c.}\hfill \\ \hfill & +ħ\delta c^{†}c+H_{\mathrm{vib}},\hfill \end{array}$$(35) ) by analysing the dynamics of the electronic states of the atoms for the parameter regime given in Equation (29(29) $$\begin{array}{cc}\hfill \mathrm{\Omega },\mathrm{\Gamma }& \gg \nu ,\kappa ,\delta ,g,\hfill \end{array}$$(29) ). Since ${\rho }_{\mathrm{A}}$ evolves according to the effective atomic master equation (30(30) $$\begin{array}{cc}\hfill {\dot{\rho }}_{\mathrm{A}}& =-\frac{\mathrm{i}\mathrm{\Omega }}{2}\left[{\sigma }_i^{-}+{\sigma }_i^{-},{\rho }_{\mathrm{A}}\right]\hfill \\ \hfill & +\sum _{i=1}^N\frac{\mathrm{\Gamma }}{2}\left(2{\sigma }_i^{-}{\rho }_{\mathrm{A}}{\sigma }_i^{+}-{\sigma }_i^{+}{\sigma }_i^{-}{\rho }_{\mathrm{A}}-{\rho }_{\mathrm{A}}{\sigma }_i^{+}{\sigma }_i^{-}\right)\hfill \end{array}$$(30) ), the time derivative of the expectation value $⟨A⟩$ of any atomic observable A equals(B1) $$\begin{array}{cc}\hfill ⟨\dot{A}⟩& =-\frac{\mathrm{i}\mathrm{\Omega }}{2}〈\left[A,{\sigma }_i^{-}+{\sigma }_i^{-}\right]〉\hfill \\ \hfill & +\sum _{i=1}^N\frac{\mathrm{\Gamma }}{2}〈\left[{\sigma }_i^{+},A\right]{\sigma }_i^{-}+{\sigma }_i^{+}\left[A,{\sigma }_i^{-}\right]〉.\hfill \end{array}$$(B1)

To obtain a closed set of rate equations, we consider the operators(B2) $$\begin{array}{c}\hfill \left({\mathrm{\Sigma }}_1^{(i)},{\mathrm{\Sigma }}_2^{(i)},{\mathrm{\Sigma }}_3^{(i)}\right)=\left({\sigma }_i^{+}{\sigma }_i^{-},{\sigma }_i^{-}+{\sigma }_i^{+},\mathrm{i}\left({\sigma }_i^{-}-{\sigma }_i^{+}\right)\right)\end{array}$$(B2)

with the commutator relations(B3) $$\begin{array}{cc}\hfill \left[{\mathrm{\Sigma }}_1^{(i)},{\mathrm{\Sigma }}_2^{(i)}\right]& =\mathrm{i}{\mathrm{\Sigma }}_3^{(i)},\left[{\mathrm{\Sigma }}_1^{(i)},{\mathrm{\Sigma }}_3^{(i)}\right]=-\mathrm{i}{\mathrm{\Sigma }}_2^{(i)},\hfill \\ \hfill \left[{\mathrm{\Sigma }}_2^{(i)},{\mathrm{\Sigma }}_3^{(i)}\right]& =-2\mathrm{i}\left(1-2{\mathrm{\Sigma }}_1^{(i)}\right)\hfill \end{array}$$(B3)

and denote their respective expectation values by $s_1$, $s_2$ and $s_3$. Combining Equations (B1(B1) $$\begin{array}{cc}\hfill ⟨\dot{A}⟩& =-\frac{\mathrm{i}\mathrm{\Omega }}{2}〈\left[A,{\sigma }_i^{-}+{\sigma }_i^{-}\right]〉\hfill \\ \hfill & +\sum _{i=1}^N\frac{\mathrm{\Gamma }}{2}〈\left[{\sigma }_i^{+},A\right]{\sigma }_i^{-}+{\sigma }_i^{+}\left[A,{\sigma }_i^{-}\right]〉.\hfill \end{array}$$(B1) )-(B3(B3) $$\begin{array}{cc}\hfill \left[{\mathrm{\Sigma }}_1^{(i)},{\mathrm{\Sigma }}_2^{(i)}\right]& =\mathrm{i}{\mathrm{\Sigma }}_3^{(i)},\left[{\mathrm{\Sigma }}_1^{(i)},{\mathrm{\Sigma }}_3^{(i)}\right]=-\mathrm{i}{\mathrm{\Sigma }}_2^{(i)},\hfill \\ \hfill \left[{\mathrm{\Sigma }}_2^{(i)},{\mathrm{\Sigma }}_3^{(i)}\right]& =-2\mathrm{i}\left(1-2{\mathrm{\Sigma }}_1^{(i)}\right)\hfill \end{array}$$(B3) ), one can then show that(B4) $$\begin{array}{cc}\hfill {\dot{s}}_1^{(i)}& =\frac{\mathrm{\Omega }}{2}{s}_3^{(i)}-\mathrm{\Gamma }{s}_1^{(i)},{\dot{s}}_2^{(i)}=-\frac{\mathrm{\Gamma }}{2}{s}_2^{(i)},\hfill \\ \hfill {\dot{s}}_3^{(i)}& =\mathrm{\Omega }\left(1-2s_1^{(i)}\right)-\frac{\mathrm{\Gamma }}{2}{s}_3^{(i)}.\hfill \end{array}$$(B4)

Consequently, the expectation values of the stationary atomic state ${\rho }_{\mathrm{A}}^{\mathrm{ss}}$ with ${\dot{s}}_1^{(i)}={\dot{s}}_2^{(i)}={\dot{s}}_3^{(i)}=0$, are given by(B5) $$\begin{array}{cc}\hfill \{s_1^{\mathrm{ss}},s_2^{\mathrm{ss}},s_3^{\mathrm{ss}}\}& =\left\{\frac{{\mathrm{\Omega }}^2}{{\mathrm{\Gamma }}^2+2{\mathrm{\Omega }}^2},0,\frac{2\mathrm{\Gamma }\mathrm{\Omega }}{{\mathrm{\Gamma }}^2+2{\mathrm{\Omega }}^2}\right\}\hfill \end{array}$$(B5)

for all particles i. This implies(B6) $$\begin{array}{cc}\hfill {}_i⟨0|{\rho }_{\mathrm{A}}^{\mathrm{ss}}{|1⟩}_i& =-\frac{\mathrm{i}\mathrm{\Gamma }\mathrm{\Omega }}{{\mathrm{\Gamma }}^2+2{\mathrm{\Omega }}^2}\hfill \end{array}$$(B6)

for all particles i, since ${}_i⟨0|{\rho }_{\mathrm{A}}^{\mathrm{ss}}{|1⟩}_i=(s_2^{\mathrm{ss}}-\mathrm{i}{s}_3^{\mathrm{ss}})/2$. Substituting the result into Equation (35(35) $$\begin{array}{cc}\hfill H_{\mathrm{eff}}& =\sum _{i=1}^Nħ\eta g{}_i⟨0|{\rho }_{\mathrm{A}}^{\mathrm{ss}}{|1⟩}_i\left(b_i+b_i^{†}\right)c+\text{H.c.}\hfill \\ \hfill & +ħ\delta c^{†}c+H_{\mathrm{vib}},\hfill \end{array}$$(35) ) yields the final the effective Hamiltonian $H_{\mathrm{eff}}$ in Equation (36(36) $$\begin{array}{cc}\hfill H_{\mathrm{eff}}& =\sum _{i=1}^N\mathrm{i}ħ\eta g_{\mathrm{eff}}\left(b_i+b_i^{†}\right)\left(c-c^{†}\right)+ħ\nu b_i^{†}{b}_i\hfill \\ \hfill & +ħ\mu {\left(b_i+b_i^{†}\right)}^3+ħ\delta c^{†}c\hfill \end{array}$$(36) ).

Relevant expectation values

In analogy to Equation (B2(B2) $$\begin{array}{c}\hfill \left({\mathrm{\Sigma }}_1^{(i)},{\mathrm{\Sigma }}_2^{(i)},{\mathrm{\Sigma }}_3^{(i)}\right)=\left({\sigma }_i^{+}{\sigma }_i^{-},{\sigma }_i^{-}+{\sigma }_i^{+},\mathrm{i}\left({\sigma }_i^{-}-{\sigma }_i^{+}\right)\right)\end{array}$$(B2) ), we now introduce the operators $B_j^{(i)}$ and $C_k$ such that(B7) $$\begin{array}{cc}\hfill \left(B_1^{(i)},B_2^{(i)},B_3^{(i)}\right)& =\left(b_i^{†}{b}_i,b_i+s{b}_i^{†},\mathrm{i}\left(b_i-b_i^{†}\right)\right),\hfill \\ \hfill \left(C_1,C_2,C_3\right)& =\left(c^{†}c,c+c^{†},\mathrm{i}\left(c-c^{†}\right)\right).\hfill \end{array}$$(B7)

Using this notation, the quantum coherence $\zeta$ and the average mean phonon number m can be written as(B8) $$\begin{array}{cc}\hfill \zeta & =\frac{1}{4N(N-1)}\sum _{i=1}^N\sum _{j\ne i}〈B_2^{(i)}{B}_2^{(j)}+B_3^{(i)}{B}_3^{(j)}〉,\hfill \\ \hfill m& =\frac{1}{N}\sum _{i=1}^N〈B_1^{(i)}〉.\hfill \end{array}$$(B8)

As we shall see below, in order to analyse the time evolution of these two variables we moreover need to consider the dynamics of the averaged expectation values(B9) $$\begin{array}{ccc}\hfill x_{\mathit{ab}}& \equiv \hfill & \hfill \frac{1}{N}\sum _{i=1}^N〈B_a^{(i)}{C}_b〉,\\ \hfill y_{\mathit{abc}}& \equiv \hfill & \hfill \frac{1}{N}\sum _{i=1}^N〈B_a^{(i)}{B}_b^{(i)}{B}_c^{(i)}〉,\\ \hfill {\tilde{y}}_{\mathit{abc}}& \equiv \hfill & \hfill \frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}〈B_a^{(i)}{B}_b^{(i)}{B}_c^{(j)}〉,\\ \hfill {\tilde{z}}_{\mathit{ab}}& \equiv \hfill & \hfill \frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}〈B_a^{(i)}{B}_b^{(j)}〉.\end{array}$$(B9)

All of the above expectation values are macroscopic quantities. They are independent of the number N of atoms in the trap in the large N limit.

Time evolution

Next we proceed as described in Section 1 and use the effective Hamiltonian $H_{\mathrm{eff}}$ in Equation (36(36) $$\begin{array}{cc}\hfill H_{\mathrm{eff}}& =\sum _{i=1}^N\mathrm{i}ħ\eta g_{\mathrm{eff}}\left(b_i+b_i^{†}\right)\left(c-c^{†}\right)+ħ\nu b_i^{†}{b}_i\hfill \\ \hfill & +ħ\mu {\left(b_i+b_i^{†}\right)}^3+ħ\delta c^{†}c\hfill \end{array}$$(36) ) together with the master equation in Equation (2(2) $$\begin{array}{c}\hfill {\dot{\rho }}_{\mathrm{I}}=-\frac{\mathrm{i}}{ħ}\left[H_{\mathrm{ideal}},{\rho }_{\mathrm{I}}\right]+\frac{\kappa }{2}\left(2c{\rho }_{\mathrm{I}}{c}^{†}-c^{†}c{\rho }_{\mathrm{I}}-{\rho }_{\mathrm{I}}{c}^{†}c\right).\end{array}$$(2) ) to derive effective rate equations for the expectation values m and $\zeta$ and the expectation values which we introduced in the previous subsection. During these calculations, we make extensive use of the commutator relations(B10) $$\begin{array}{cc}\hfill \left[B_1^{(i)},B_2^{(i)}\right]& =\mathrm{i}{B}_3^{(i)},\left[B_1^{(i)},B_3^{(i)}\right]=-\mathrm{i}{B}_2^{(i)},\hfill \\ \hfill \left[B_2^{(i)},B_3^{(i)}\right]& =-2\mathrm{i},[C_1,C_2]=\mathrm{i}{C}_3,\hfill \\ \hfill [C_1,C_3]& =-\mathrm{i}{C}_2,[C_2,C_3]=-2\mathrm{i}\hfill \end{array}$$(B10)

which we summarise here for convenience. Moreover, using Equation (B7(B7) $$\begin{array}{cc}\hfill \left(B_1^{(i)},B_2^{(i)},B_3^{(i)}\right)& =\left(b_i^{†}{b}_i,b_i+s{b}_i^{†},\mathrm{i}\left(b_i-b_i^{†}\right)\right),\hfill \\ \hfill \left(C_1,C_2,C_3\right)& =\left(c^{†}c,c+c^{†},\mathrm{i}\left(c-c^{†}\right)\right).\hfill \end{array}$$(B7) ), Equation (15(15) $$\begin{array}{cc}\hfill H_{\mathrm{I}}& =\sum _{i=1}^N\frac{ħ\mathrm{\Omega }}{2}{\sigma }_i^{-}+ħ\eta g(b_i+b_i^{†})c{\sigma }_i^{+}+\text{H.c.}\hfill \\ \hfill & +ħ\delta c^{†}c+H_{\mathrm{vib}}.\hfill \end{array}$$(15) ) simplifies to(B11) $$\begin{array}{cc}\hfill H_{\mathrm{eff}}& =\sum _{i=1}^Nħ\eta g_{\mathrm{eff}}{B}_2^{(i)}{C}_3+ħ\nu B_1^{(i)}+ħ\mu B_2^{(i)3}\hfill \\ \hfill & +ħ\delta C_1.\hfill \end{array}$$(B11)

In the following, we combine Equation (A1(A1) $$\begin{array}{cc}\hfill ⟨\dot{A}⟩& =-\frac{\mathrm{i}}{ħ}〈\left[A,H_{\mathrm{ideal}}\right]〉+\frac{\kappa }{2}〈[c^{†},A]c+c^{†}[A,c]〉.\hfill \end{array}$$(A1) ) (but with $H_{\mathrm{ideal}}$ replaced by $H_{\mathrm{eff}}$) with Equations (B10(B10) $$\begin{array}{cc}\hfill \left[B_1^{(i)},B_2^{(i)}\right]& =\mathrm{i}{B}_3^{(i)},\left[B_1^{(i)},B_3^{(i)}\right]=-\mathrm{i}{B}_2^{(i)},\hfill \\ \hfill \left[B_2^{(i)},B_3^{(i)}\right]& =-2\mathrm{i},[C_1,C_2]=\mathrm{i}{C}_3,\hfill \\ \hfill [C_1,C_3]& =-\mathrm{i}{C}_2,[C_2,C_3]=-2\mathrm{i}\hfill \end{array}$$(B10) ) and (B11(B11) $$\begin{array}{cc}\hfill H_{\mathrm{eff}}& =\sum _{i=1}^Nħ\eta g_{\mathrm{eff}}{B}_2^{(i)}{C}_3+ħ\nu B_1^{(i)}+ħ\mu B_2^{(i)3}\hfill \\ \hfill & +ħ\delta C_1.\hfill \end{array}$$(B11) ) to obtain the time derivatives in Equation (38(38) $$\begin{array}{cc}\hfill \dot{m}& =\eta g_{\mathrm{eff}}{x}_{33}+3\mu y_{232},\hfill \\ \hfill \dot{\zeta }& =\eta g_{\mathrm{eff}}{x}_{33}+3\mu {\tilde{y}}_{223}\hfill \end{array}$$(38) ) which describe the relatively slow dynamics of m and $\zeta$.

Given the parameter regime in Equation (39(39) $$\begin{array}{cc}\hfill N\eta g_{\mathrm{eff}},\nu ,\kappa ,\delta & \gg \eta g_{\mathrm{eff}},\mu .\hfill \end{array}$$(39) ), the expectation values $x_{\mathit{ab}}$ which we defined in Equation (B9(B9) $$\begin{array}{ccc}\hfill x_{\mathit{ab}}& \equiv \hfill & \hfill \frac{1}{N}\sum _{i=1}^N〈B_a^{(i)}{C}_b〉,\\ \hfill y_{\mathit{abc}}& \equiv \hfill & \hfill \frac{1}{N}\sum _{i=1}^N〈B_a^{(i)}{B}_b^{(i)}{B}_c^{(i)}〉,\\ \hfill {\tilde{y}}_{\mathit{abc}}& \equiv \hfill & \hfill \frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}〈B_a^{(i)}{B}_b^{(i)}{B}_c^{(j)}〉,\\ \hfill {\tilde{z}}_{\mathit{ab}}& \equiv \hfill & \hfill \frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}〈B_a^{(i)}{B}_b^{(j)}〉.\end{array}$$(B9) ) evolve on a much faster time scale than m and $\zeta$. Taking this into account and neglecting terms which scale as $\eta g_{\mathrm{eff}}$ and $\mu$, we find that the single-particle expectation values $x_{\mathit{ab}}$ evolve to a very good approximation according to the differential equations(B12) $$\begin{array}{cc}\hfill {\dot{x}}_{22}& =-2N\eta g_{\mathrm{eff}}{\tilde{z}}_{22}-\nu x_{32}-\delta x_{23}-\frac{\kappa }{2}{x}_{22},\hfill \\ \hfill {\dot{x}}_{23}& =-\nu x_{33}+\delta x_{22}-\frac{\kappa }{2}{x}_{23},\hfill \\ \hfill {\dot{x}}_{32}& =-2N\eta g_{\mathrm{eff}}{\tilde{z}}_{23}+\nu x_{22}-\delta x_{33}-\frac{\kappa }{2}{x}_{32},\hfill \\ \hfill {\dot{x}}_{33}& =\nu x_{23}+\delta x_{32}-\frac{\kappa }{2}{x}_{33}.\hfill \end{array}$$(B12)

Moreover, the time derivative of the two-particle quantum coherences ${\tilde{z}}_{\mathit{ab}}$ are given by(B13) $$\begin{array}{cc}\hfill {\dot{\tilde{z}}}_{22}& =-2\nu {\tilde{z}}_{23},\hfill \\ \hfill {\dot{\tilde{z}}}_{23}={\dot{\tilde{z}}}_{32}& =\nu \left({\tilde{z}}_{22}-{\tilde{z}}_{33}\right),\hfill \\ \hfill {\dot{\tilde{z}}}_{33}& =2\nu {\tilde{z}}_{23},\hfill \end{array}$$(B13)

if we neglect again terms which scale as $\eta g_{\mathrm{eff}}$ or $\mu$. Equation (B13(B13) $$\begin{array}{cc}\hfill {\dot{\tilde{z}}}_{22}& =-2\nu {\tilde{z}}_{23},\hfill \\ \hfill {\dot{\tilde{z}}}_{23}={\dot{\tilde{z}}}_{32}& =\nu \left({\tilde{z}}_{22}-{\tilde{z}}_{33}\right),\hfill \\ \hfill {\dot{\tilde{z}}}_{33}& =2\nu {\tilde{z}}_{23},\hfill \end{array}$$(B13) ) describes relatively fast oscillations on the time scale given by the phonon frequency $\nu$. As we shall see in the next subsection of this appendix, on the relatively slow time scale we can therefore approximate these coherences, which depend on the slowly evolving variable $\zeta$, by their average values.

To obtain a closed set of effective cooling equations, we also need to have a closer look at the y and the $\tilde{y}$ expectation values in Equation (38(38) $$\begin{array}{cc}\hfill \dot{m}& =\eta g_{\mathrm{eff}}{x}_{33}+3\mu y_{232},\hfill \\ \hfill \dot{\zeta }& =\eta g_{\mathrm{eff}}{x}_{33}+3\mu {\tilde{y}}_{223}\hfill \end{array}$$(38) ). Doing so, we find, for example, that(B14) $$\begin{array}{cc}\hfill {\dot{y}}_{222}& =-3\nu y_{232},\hfill \end{array}$$(B14)

while(B15) $$\begin{array}{cc}\hfill {\dot{\tilde{y}}}_{222}& =-\nu \left({\tilde{y}}_{223}+2{\tilde{y}}_{232}\right),\hfill \\ \hfill {\dot{\tilde{y}}}_{233}& =\nu \left({\tilde{y}}_{223}+{\tilde{y}}_{232}-{\tilde{y}}_{333}\right),\hfill \\ \hfill {\dot{\tilde{y}}}_{332}& =\nu \left(2{\tilde{y}}_{232}-{\tilde{y}}_{333}\right).\hfill \end{array}$$(B15)

These expectation values all oscillate essentially rapidly around zero.

Adiabatic elimination

To finally obtain an effective cooling equation for the time evolution of the slowly evolving phonon expectation values m and $\zeta$, we now determine the x, y, $\tilde{y}$ and the $\tilde{z}$-expectation values in Equation (38(38) $$\begin{array}{cc}\hfill \dot{m}& =\eta g_{\mathrm{eff}}{x}_{33}+3\mu y_{232},\hfill \\ \hfill \dot{\zeta }& =\eta g_{\mathrm{eff}}{x}_{33}+3\mu {\tilde{y}}_{223}\hfill \end{array}$$(38) ) via an adiabatic elimination. Considering the parameter regime in Equation (39(39) $$\begin{array}{cc}\hfill N\eta g_{\mathrm{eff}},\nu ,\kappa ,\delta & \gg \eta g_{\mathrm{eff}},\mu .\hfill \end{array}$$(39) ) and setting the time derivatives in Equation (B13(B13) $$\begin{array}{cc}\hfill {\dot{\tilde{z}}}_{22}& =-2\nu {\tilde{z}}_{23},\hfill \\ \hfill {\dot{\tilde{z}}}_{23}={\dot{\tilde{z}}}_{32}& =\nu \left({\tilde{z}}_{22}-{\tilde{z}}_{33}\right),\hfill \\ \hfill {\dot{\tilde{z}}}_{33}& =2\nu {\tilde{z}}_{23},\hfill \end{array}$$(B13) ) equal to zero, we find for example that(B16) $$\begin{array}{cc}\hfill {\tilde{z}}_{23}={\tilde{z}}_{32}& =0,\hfill \\ \hfill {\tilde{z}}_{22}& ={\tilde{z}}_{33}.\hfill \end{array}$$(B16)

Comparing $\zeta$ in Equation (B8(B8) $$\begin{array}{cc}\hfill \zeta & =\frac{1}{4N(N-1)}\sum _{i=1}^N\sum _{j\ne i}〈B_2^{(i)}{B}_2^{(j)}+B_3^{(i)}{B}_3^{(j)}〉,\hfill \\ \hfill m& =\frac{1}{N}\sum _{i=1}^N〈B_1^{(i)}〉.\hfill \end{array}$$(B8) ) with the definition of ${\tilde{z}}_{\mathit{ab}}$ in Equation (B9(B9) $$\begin{array}{ccc}\hfill x_{\mathit{ab}}& \equiv \hfill & \hfill \frac{1}{N}\sum _{i=1}^N〈B_a^{(i)}{C}_b〉,\\ \hfill y_{\mathit{abc}}& \equiv \hfill & \hfill \frac{1}{N}\sum _{i=1}^N〈B_a^{(i)}{B}_b^{(i)}{B}_c^{(i)}〉,\\ \hfill {\tilde{y}}_{\mathit{abc}}& \equiv \hfill & \hfill \frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}〈B_a^{(i)}{B}_b^{(i)}{B}_c^{(j)}〉,\\ \hfill {\tilde{z}}_{\mathit{ab}}& \equiv \hfill & \hfill \frac{1}{N(N-1)}\sum _{i=1}^N\sum _{j\ne i}〈B_a^{(i)}{B}_b^{(j)}〉.\end{array}$$(B9) ), we moreover find that $\zeta =({\tilde{z}}_{22}+{\tilde{z}}_{33})/4$ which implies(B17) $$\begin{array}{cc}\hfill {\tilde{z}}_{22}={\tilde{z}}_{33}& =2\zeta .\hfill \end{array}$$(B17)

Taking this into account when adiabatically eliminating the $x_{\mathit{ab}}$ variables by setting the time derivatives in Equation (B12(B12) $$\begin{array}{cc}\hfill {\dot{x}}_{22}& =-2N\eta g_{\mathrm{eff}}{\tilde{z}}_{22}-\nu x_{32}-\delta x_{23}-\frac{\kappa }{2}{x}_{22},\hfill \\ \hfill {\dot{x}}_{23}& =-\nu x_{33}+\delta x_{22}-\frac{\kappa }{2}{x}_{23},\hfill \\ \hfill {\dot{x}}_{32}& =-2N\eta g_{\mathrm{eff}}{\tilde{z}}_{23}+\nu x_{22}-\delta x_{33}-\frac{\kappa }{2}{x}_{32},\hfill \\ \hfill {\dot{x}}_{33}& =\nu x_{23}+\delta x_{32}-\frac{\kappa }{2}{x}_{33}.\hfill \end{array}$$(B12) ) equal to zero yields(B18) $$\begin{array}{cc}\hfill x_{22}& =-\frac{8N\eta g_{\mathrm{eff}}\kappa \left[{\kappa }^2+4\left({\delta }^2+{\nu }^2\right)\right]}{{\kappa }^4+16{\left({\delta }^2-{\nu }^2\right)}^2+8{\kappa }^2\left({\delta }^2+{\nu }^2\right)}\zeta ,\hfill \\ \hfill x_{23}& =-\frac{16N\eta g_{\mathrm{eff}}\delta \left[{\kappa }^2+4\left({\delta }^2-{\nu }^2\right)\right]}{{\kappa }^4+16{\left({\delta }^2-{\nu }^2\right)}^2+8{\kappa }^2\left({\delta }^2+{\nu }^2\right)}\zeta ,\hfill \\ \hfill x_{32}& =-\frac{16N\eta g_{\mathrm{eff}}\nu \left[{\kappa }^2-4\left({\delta }^2-4{\nu }^2\right)\right]}{{\kappa }^4+16{\left({\delta }^2-{\nu }^2\right)}^2+8{\kappa }^2\left({\delta }^2+{\nu }^2\right)}\zeta ,\hfill \\ \hfill x_{33}& =-\frac{64N\eta g_{\mathrm{eff}}\kappa \delta \nu }{{\kappa }^4+16{\left({\delta }^2-{\nu }^2\right)}^2+8{\kappa }^2\left({\delta }^2+{\nu }^2\right)}\zeta .\hfill \end{array}$$(B18)

This result can be checked relatively easily by substituting Equations (B17(B17) $$\begin{array}{cc}\hfill {\tilde{z}}_{22}={\tilde{z}}_{33}& =2\zeta .\hfill \end{array}$$(B17) ) and (B18(B18) $$\begin{array}{cc}\hfill x_{22}& =-\frac{8N\eta g_{\mathrm{eff}}\kappa \left[{\kappa }^2+4\left({\delta }^2+{\nu }^2\right)\right]}{{\kappa }^4+16{\left({\delta }^2-{\nu }^2\right)}^2+8{\kappa }^2\left({\delta }^2+{\nu }^2\right)}\zeta ,\hfill \\ \hfill x_{23}& =-\frac{16N\eta g_{\mathrm{eff}}\delta \left[{\kappa }^2+4\left({\delta }^2-{\nu }^2\right)\right]}{{\kappa }^4+16{\left({\delta }^2-{\nu }^2\right)}^2+8{\kappa }^2\left({\delta }^2+{\nu }^2\right)}\zeta ,\hfill \\ \hfill x_{32}& =-\frac{16N\eta g_{\mathrm{eff}}\nu \left[{\kappa }^2-4\left({\delta }^2-4{\nu }^2\right)\right]}{{\kappa }^4+16{\left({\delta }^2-{\nu }^2\right)}^2+8{\kappa }^2\left({\delta }^2+{\nu }^2\right)}\zeta ,\hfill \\ \hfill x_{33}& =-\frac{64N\eta g_{\mathrm{eff}}\kappa \delta \nu }{{\kappa }^4+16{\left({\delta }^2-{\nu }^2\right)}^2+8{\kappa }^2\left({\delta }^2+{\nu }^2\right)}\zeta .\hfill \end{array}$$(B18) ) back into Equation (B12(B12) $$\begin{array}{cc}\hfill {\dot{x}}_{22}& =-2N\eta g_{\mathrm{eff}}{\tilde{z}}_{22}-\nu x_{32}-\delta x_{23}-\frac{\kappa }{2}{x}_{22},\hfill \\ \hfill {\dot{x}}_{23}& =-\nu x_{33}+\delta x_{22}-\frac{\kappa }{2}{x}_{23},\hfill \\ \hfill {\dot{x}}_{32}& =-2N\eta g_{\mathrm{eff}}{\tilde{z}}_{23}+\nu x_{22}-\delta x_{33}-\frac{\kappa }{2}{x}_{32},\hfill \\ \hfill {\dot{x}}_{33}& =\nu x_{23}+\delta x_{32}-\frac{\kappa }{2}{x}_{33}.\hfill \end{array}$$(B12) ). Finally, setting the right hand sides of Equations (B14(B14) $$\begin{array}{cc}\hfill {\dot{y}}_{222}& =-3\nu y_{232},\hfill \end{array}$$(B14) ) and (B15(B15) $$\begin{array}{cc}\hfill {\dot{\tilde{y}}}_{222}& =-\nu \left({\tilde{y}}_{223}+2{\tilde{y}}_{232}\right),\hfill \\ \hfill {\dot{\tilde{y}}}_{233}& =\nu \left({\tilde{y}}_{223}+{\tilde{y}}_{232}-{\tilde{y}}_{333}\right),\hfill \\ \hfill {\dot{\tilde{y}}}_{332}& =\nu \left(2{\tilde{y}}_{232}-{\tilde{y}}_{333}\right).\hfill \end{array}$$(B15) ) equal to zero, yields(B19) $$\begin{array}{cc}\hfill y_{232}={\tilde{y}}_{223}& =0\hfill \end{array}$$(B19)

up to terms which scale as $\mu /\nu$. Combining the results of this subsection with Equation (38(38) $$\begin{array}{cc}\hfill \dot{m}& =\eta g_{\mathrm{eff}}{x}_{33}+3\mu y_{232},\hfill \\ \hfill \dot{\zeta }& =\eta g_{\mathrm{eff}}{x}_{33}+3\mu {\tilde{y}}_{223}\hfill \end{array}$$(38) ) eventually yields the effective cooling dynamics in Equation (40(40) $$\begin{array}{cc}\hfill \dot{m}=\dot{\zeta }& =-A_N\zeta \hfill \end{array}$$(40) ).