A quantum optomechanical interface beyond the resolved sideband limit

Quantum interfaces are anticipated to form a crucial component of future quantum information networks due to their ability to transfer quantum information between 'flying' carriers—photons—and 'stationary' quantum media [1, 2]. Numerous physical settings lend themselves to creating such interfaces, including both warm and cold neutral atoms [3, 4], quantum dots [5, 6], nitrogen-vacancy (and other) colour centres [7, 8] and trapped ions [9, 10].

Cavity optomechanical systems are a promising platform for novel quantum interfaces because they may be used to couple light to the 'stationary' media listed above, and because existing mechanical oscillators have extremely long decay times (∼30 s, e.g. [11]). In these systems a mechanical degree of freedom with (angular) frequency ${\omega }_{{\rm{M}}}$ is parametrically coupled to an optical cavity with linewidth κ. These devices have entered the quantum regime for the first time in recent years (e.g. [12–15]) and they are set to find technological and research applications in a variety of capacities. For instance, a mechanical oscillator which is subject to radiation pressure arising from microwave and optical modes simultaneously may be used to couple superconducting circuits to optics [16–21]; the capability to do this at the quantum level would permit construction of quantum information networks that utilise the advantages of both low-loss optical transmission and the exquisite quantum control of superconducting circuits. There are also a wealth of proposals for employing cavity optomechanical coupling to prepare massive mechanical oscillators in non-classical states, which could allow sensitive probing of the quantum-to-classical transition [12] or quantum-enhanced metrology (e.g. mass sensing [22], accelerometry [23]).

In many of these applications it is expedient to employ a low-frequency mechanical oscillator (e.g. [24]); for instance, devices with material-limited damping rates Γ obey the scaling ${\rm{\Gamma }}\propto {\omega }_{{\rm{M}}}^{2}$ (due to the Akhiezer effect [25]). Furthermore, large optomechanical coupling rates require small cavity volumes, and the optical linewidth typically scales inversely with its length, so that κ is often unavoidably larger than ${\omega }_{{\rm{M}}}$ in micro- and nano-optomechanical devices. This high-bandwidth operating regime, $\kappa \gg {\omega }_{{\rm{M}}}$, is known as the unresolved sideband limit or bad-cavity limit, despite the fact that the cavities involved often have large Q factors.

Current technologies make it possible to cool a mechanical oscillator to near its ground state in the unresolved sideband limit using optical measurement-based feedback cooling (e.g. [26]), or to perform steady-state cooling using hybrid quantum systems [27–29], dissipative optomechanics [30], optomechanically induced transparency [31], and related schemes (e.g. [32, 33]). Very rapid ground state cooling of the mechanical oscillator is also enabled by varying the optical drive's amplitude and phase dynamically [34–38]. Alternatively, one may exploit the large cavity bandwidth in the bad-cavity limit to engineer quantum non-demolition (QND) interactions between short optical pulses and a mechanical oscillator (see (5)) [34]. QND interactions are the basis of proposed quantum 'upload' interfaces, which write the state of one input mode onto one output mode [39–41]. Further applications of QND interactions include generation of mechanical geometric phases [35], and cooling by measurement and state tomography [36].

Despite this—and the success of optomechanical interfaces in the 'resolved-sideband' ($\kappa \ll {\omega }_{{\rm{M}}}$), regime (e.g. [14, 15, 42–45])—there are currently no quantum optical-to-mechanical interfaces which permit a state swap in the unresolved sideband limit. Teleportation schemes [46] and the aforementioned quantum 'upload' interfaces [39, 40], which transfer an input state from one degree of freedom to another, are incapable of directly exchanging two input states.

Here we present a bidirectional optomechanical interface which can perform a quantum state swap and other operations in the unresolved sideband limit. Our interface is completely deterministic and does not rely on measurement or any form of conditioning or post-selection. Furthermore, the proposed protocol only requires classical open-loop control, which may be performed with very low added noise (i.e. below the back-action noise) [47].

We predict that the interface is capable of performing a number of important operations, including mechanical ground state cooling, squeezing, and exploiting non-Gaussian optical inputs to create Wigner-negative mechanical states, including low-amplitude Schrödinger cat ('kitten') states.

We now show that QND interactions may be performed by pulsed optomechanics, following [34]. Consider a mechanical oscillator with annihilation operator $\hat{b}$ and angular frequency ${\omega }_{{\rm{M}}}$ parametrically coupled to an optical cavity. As the (dimensionless) position of the oscillator ${X}_{{\rm{M}}}={\hat{b}}^{\dagger }+\hat{b}$ (momentum ${P}_{{\rm{M}}}={\rm{i}}({\hat{b}}^{\dagger }-\hat{b})$) fluctuates the resonance frequency of the cavity also shifts. Thus the optomechanical (dispersive) coupling Hamiltonian takes the form [24, 48]

$$\begin{eqnarray*}&&\hat{H}={\hslash }{g}_{0}{X}_{{\rm{M}}}{\hat{a}}^{\dagger }\hat{a},\end{eqnarray*}$$

where $\hat{a}$ is the cavity annihilation operator and the optomechanical coupling rate is g₀. This is understood as the frequency shift of the cavity per zero-point displacement of the oscillator. For a detailed discussion of the optomechanical Hamiltonian the reader is referred to [24].

We assume that ${g}_{0}\ll {\omega }_{{\rm{M}}}$, as is usually the case (e.g. review of parameters given by [49]). Then if a bright pump beam is applied to the cavity the interaction between the cavity field and mechanical fluctuations is well described by the effective linearised Hamiltonian [34]

$$\begin{eqnarray*}&&\hat{H}={\hslash }{g}_{0}| \alpha (t)| {X}_{{\rm{M}}}{\bar{X}}_{{\rm{L}}}(t)+{\hslash }{g}_{0}{X}_{{\rm{M}}}{| \alpha (t)| }^{2}.\end{eqnarray*}$$

Here, $\alpha (t)$ is the mean intracavity amplitude induced by the pump, and the optical fluctuation operator ${\bar{X}}_{{\rm{L}}}(t)={\hat{a}}^{\dagger }+\hat{a}-2\alpha (t)$ describes amplitude fluctuations of the field (this expression is valid for real, positive α; for the general case see appendix A). These quantities are written with explicit time dependence to show that they depend on the pulse used to illuminate the cavity. The phase fluctuations of the optical field ${\bar{P}}_{{\rm{L}}}(t)$ do not appear in the interaction Hamiltonian. Both fluctuation operators obey $\langle {\bar{X}}_{{\rm{L}}}(t)\rangle =\langle {\bar{P}}_{{\rm{L}}}(t)\rangle =0$ for all times.

We now consider applying pulses of light to the cavity. These are described by a mean amplitude ${\alpha }_{\mathrm{in}}(t)$, normalised such that the average number of photons in the input pulse envelope is

$$\begin{eqnarray*}&&N={\int }_{-\infty }^{\infty }{\rm{d}}t\ \ {| {\alpha }_{\mathrm{in}}(t)| }^{2}\gg 1\end{eqnarray*}$$

and the instantaneous amplitude and phase noise operators, ${\bar{X}}_{{\rm{L}},\mathrm{in}}(t)$ and ${\bar{P}}_{{\rm{L}},\mathrm{in}}(t)$. The latter are defined by analogy with the intracavity field. The optical field exiting the cavity is connected to the intracavity mode and input field through the input–output relation. This takes the form

$$\begin{eqnarray}&&{\bar{X}}_{{\rm{L}},\mathrm{out}}(t)=\sqrt{\kappa }{\bar{X}}_{{\rm{L}}}(t)-{\bar{X}}_{{\rm{L}},\mathrm{in}}(t),\end{eqnarray} \tag{ 1 }$$

with a similar expression holding for the phase fluctuations [50].

If the bandwidth of the pulse is large compared to both the mechanical frequency and the thermal heating rate ${\rm{\Gamma }}\times ({\bar{n}}_{{\rm{B}}}+1/2)$ we may neglect the damped dynamics of the oscillator over the duration of the interaction; furthermore, the large cavity linewidth (compared to the pulse bandwidth) allows the optical fields to be adiabatically eliminated. Together, these approximations yield effective Langevin equations of the form

$$\begin{eqnarray}&&{\dot{X}}_{{\rm{M}}}=0,\quad \end{eqnarray} \tag{ 2a }$$

$$\begin{eqnarray}&&{\dot{P}}_{{\rm{M}}}=-\displaystyle \frac{8{g}_{0}| {\alpha }_{\mathrm{in}}(t)| }{\kappa }({\bar{X}}_{{\rm{L}},\mathrm{in}}(t)+| {\alpha }_{\mathrm{in}}(t)| ),\phantom{\rule{-6.78em}{0ex}}\end{eqnarray} \tag{ 2b }$$

$$\begin{eqnarray}&&{\bar{X}}_{{\rm{L}},\mathrm{out}}(t)={\bar{X}}_{{\rm{L}},\mathrm{in}}(t),\end{eqnarray} \tag{ 2c }$$

$$\begin{eqnarray}&&{\bar{P}}_{{\rm{L}},\mathrm{out}}(t)={\bar{P}}_{{\rm{L}},\mathrm{in}}(t)-\displaystyle \frac{8{g}_{0}| {\alpha }_{\mathrm{in}}(t)| }{\kappa }{X}_{{\rm{M}}},\phantom{\rule{-4em}{0ex}}\end{eqnarray} \tag{ 2d }$$

in which we have used (1) to obtain (2c) and (2d).

Integrating (2a) and (2b) over the duration of the pulse and employing the fact that the pulse envelope changes rapidly compared to ${\omega }_{{\rm{M}}}^{-1}$ allows us to obtain the following transformation between the initial and final mechanical quadratures (final quadratures are primed)

$$\begin{eqnarray}{X}_{{\rm{M}}}^{\prime } & = & {X}_{{\rm{M}}},\\ {P}_{{\rm{M}}}^{\prime } & = & {P}_{{\rm{M}}}+\chi \sqrt{N}+\chi {\int }_{-\infty }^{\infty }{\rm{d}}t\ \ \displaystyle \frac{| {\alpha }_{\mathrm{in}}(t)| }{\sqrt{N}}{\bar{X}}_{{\rm{L}},\mathrm{in}}(t),\end{eqnarray} \tag{ 3 }$$

where the interaction strength is $\chi =-8{g}_{0}\sqrt{N}/\kappa$. We see that the weighted optical amplitude fluctuations have been transferred onto the mechanical momentum ${P}_{{\rm{M}}}^{\prime }$, along with a classical momentum kick $\chi \sqrt{N}$. It is possible to subtract off the classical kick by performing a displacement on the mechanical oscillator immediately after the interaction; this can be done using open-loop control with very little added noise (e.g. [47]). We will therefore neglect this kick hereafter.

Motivated by the form of (3), we factorise the pulse envelope into the form ${\alpha }_{\mathrm{in}}(t)=\sqrt{N}f(t);$ thus the integral of ${| f(t)| }^{2}$ over the whole pulse is unity. Multiplying (2c) and (2d) by $| f(t)|$ and integrating, we find

$$\begin{eqnarray}&&{X}_{{\rm{L}}}^{\prime }={X}_{{\rm{L}}},\ {P}_{{\rm{L}}}^{\prime }={P}_{{\rm{L}}}+\chi {X}_{{\rm{M}}},\end{eqnarray} \tag{ 4 }$$

in which the collective pulse quadratures X_L and P_L ($[{X}_{{\rm{L}}},{P}_{{\rm{L}}}]=[{X}_{{\rm{M}}},{P}_{{\rm{M}}}]=2{\rm{i}}$) have been defined by

$$\begin{eqnarray*}{X}_{{\rm{L}}} & = & {\int }_{-\infty }^{\infty }{\rm{d}}t\ \ | f(t)| {\bar{X}}_{{\rm{L}},\mathrm{in}}(t),\\ {P}_{{\rm{L}}} & = & {\int }_{-\infty }^{\infty }{\rm{d}}t\ \ | f(t)| {\bar{P}}_{{\rm{L}},\mathrm{in}}(t).\end{eqnarray*}$$

Note that these quadratures are zero-mean by construction. Similar relations hold for ${X}_{{\rm{L}}}^{\prime }$ and ${P}_{{\rm{L}}}^{\prime }$, with the input fields replaced by ${\bar{X}}_{{\rm{L}},\mathrm{out}}(t)$ and ${\bar{P}}_{{\rm{L}},\mathrm{out}}(t)$ respectively.

The evolution described by (3) (including a mechanical displacement) and (4) is thus a standard QND interaction between the mechanical oscillator and the collective quadratures of the optical pulse;

$$\begin{eqnarray}&&{X}_{{\rm{M}}}^{\prime }={X}_{{\rm{M}}},{P}_{{\rm{M}}}^{\prime }={P}_{{\rm{M}}}+\chi {X}_{{\rm{L}}},\end{eqnarray} \tag{ 5a }$$

$$\begin{eqnarray}&&{X}_{{\rm{L}}}^{\prime }={X}_{{\rm{L}}},{P}_{{\rm{L}}}^{\prime }={P}_{{\rm{L}}}+\chi {X}_{{\rm{M}}}.\end{eqnarray} \tag{ 5b }$$

Each oscillator acquires information about one quadrature of the other during this interaction.

A QND interaction correlates each oscillator with the other, but only one quadrature at a time; evidently some extra ingredients will be required to perform a complete state swap. We will first show how to use a sequence of QND interactions and local operations to realise a state swap under ideal conditions. Secondly, we will generalise our protocol to allow both state transfer and squeezing. In our approach the pulse enters and leaves the optomechanical cavity multiple times, with local operations on the pulse and mechanical resonator occurring between the optomechanical interactions. Figure 1 summarises the necessary steps, as presented below, graphically.

Zoom In Zoom Out Reset image size

This approach is in the spirit of the proposal of [39, 40]; however, our protocol replaces measurement and feedback with a third QND interaction to achieve a two-way interface.

Firstly, consider the case in which $\chi =-1$ and there are no optical or mechanical decoherence mechanisms. From (5) it is clear that a single QND pulse correlates ${P}_{{\rm{M}}}^{\prime }$ (${P}_{{\rm{L}}}^{\prime }$) with X_L (X_M). In order to build the necessary correlations between the other pair of quadratures we must exchange ${X}_{{\rm{M}}}^{\prime }$ (${X}_{{\rm{L}}}^{\prime }$) and ${P}_{{\rm{M}}}^{\prime }$ (${P}_{{\rm{L}}}^{\prime }$). This may be achieved by local $\pi /2$ rotations on each mode.

For definiteness, we can consider the local rotations to be generated in the following manner. Delay the pulse for a quarter of a mechanical period before re-injecting it into the cavity to achieve a mechanical rotation of $\pi /2$. The optical amplitude and phase noise rotations may be generated by displacing the optical state in phase space, as demonstrated in figure 1 (panel (ii)), e.g. by interfering the pulse with a very bright, temporally matched pulse on a highly asymmetric beamsplitter (see figure 1, panel (iii)).

If we apply these rotations and then allow the pulse and oscillator to interact for a second time the quadratures have been transformed according to

$$\begin{eqnarray}&&{X}_{{\rm{M}}}^{\prime }={P}_{{\rm{M}}}-{X}_{{\rm{L}}},\;{P}_{{\rm{M}}}^{\prime }=-{P}_{{\rm{L}}},\end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}&&{X}_{{\rm{L}}}^{\prime }={P}_{{\rm{L}}}-{X}_{{\rm{M}}},\;{P}_{{\rm{L}}}^{\prime }=-{P}_{{\rm{M}}}.\end{eqnarray} \tag{ 6b }$$

The momentum quadratures have been entirely swapped, as desired, because the contribution of the initial mechanical position to ${P}_{{\rm{M}}}^{\prime }$ due to free evolution has been coherently cancelled by the contribution of ${X}_{{\rm{L}}}^{\prime }$ (which depends on the initial position) from the QND interaction. However, it is clear that ${X}_{{\rm{M}}}^{\prime }$ (${X}_{{\rm{L}}}^{\prime }$) still retains some 'memory' of the initial momentum P_M (P_L) of the mechanical (optical) oscillator; this needs to be erased for the swap to be completed. Fortunately, this can be achieved by a second pair of $\pi /2$ rotations and a third QND interaction, viz

$$\begin{eqnarray}&&{X}_{{\rm{M}}}^{\prime }=-{P}_{{\rm{L}}},{P}_{{\rm{M}}}^{\prime }={X}_{{\rm{L}}},\end{eqnarray} \tag{ 7a }$$

$$\begin{eqnarray}&&{X}_{{\rm{L}}}^{\prime }=-{P}_{{\rm{M}}},{P}_{{\rm{L}}}^{\prime }={X}_{{\rm{M}}}.\end{eqnarray} \tag{ 7b }$$

This is a perfect state swap, up to local rotations (which we shall not consider), because the mechanical system now carries only fluctuations from the initial optical state, and vice versa.

The fundamental limits on the interface's performance are imposed by optical loss and thermomechanical noise.

Optical loss—such as coupling inefficiencies, scattering and absorption—is readily treated by coupling the circulating pulse to a vacuum bath on a beamsplitter with efficiency ${\eta }_{{\rm{L}}}$ (unity for perfect operation). The resulting nonhomogeneous evolution is

$$\begin{eqnarray*}&&{X}_{{\rm{L}}}\to \sqrt{{\eta }_{{\rm{L}}}}{X}_{{\rm{L}}}+\sqrt{1-{\eta }_{{\rm{L}}}}{X}_{{\rm{V}}},\end{eqnarray*}$$

where X_V is vacuum noise. An analogous relation holds for P_L. The properties of the vacuum mode are $\langle {X}_{{\rm{V}}}^{2}\rangle =\langle {P}_{{\rm{V}}}^{2}\rangle =1$ and ${\langle {X}_{{\rm{V}}}{P}_{{\rm{V}}}\rangle }_{{\rm{S}}}=0$. Here we have introduced the symmetrised expectation value, ${\langle {X}_{{\rm{V}}}{P}_{{\rm{V}}}\rangle }_{{\rm{S}}}=\tfrac{1}{2}\langle {X}_{{\rm{V}}}{P}_{{\rm{V}}}+{P}_{{\rm{V}}}{X}_{{\rm{V}}}\rangle ={\mathfrak{R}}\{\langle {X}_{{\rm{V}}}{P}_{{\rm{V}}}\rangle \}$.

Mechanical decoherence is induced by two intimately related mechanisms: loss of energy from the oscillator into a heat bath, and random excitation of the oscillator by the bath. These processes may be described by the Langevin equations for the damped evolution (see [51])

$$\begin{eqnarray}&&{\dot{X}}_{{\rm{M}}}=+{\omega }_{{\rm{M}}}{P}_{{\rm{M}}},\end{eqnarray} \tag{ 8a }$$

$$\begin{eqnarray}&&{\dot{P}}_{{\rm{M}}}=-{\omega }_{{\rm{M}}}{X}_{{\rm{M}}}+\sqrt{2{\rm{\Gamma }}}\hat{\xi }-{\rm{\Gamma }}{P}_{{\rm{M}}}.\end{eqnarray} \tag{ 8b }$$

Momentum-dependent damping occurs at a rate Γ while excitations enter via the random force $\hat{\xi }$ at a rate of ${\rm{\Gamma }}\times ({\bar{n}}_{{\rm{B}}}+1/2)$, where ${\bar{n}}_{{\rm{B}}}$ is the equilibrium occupancy of the oscillator, related to the bath temperature T through the Bose–Einstein distribution (hence ${\bar{n}}_{{\rm{B}}}\approx {k}_{{\rm{B}}}T/{\hslash }{\omega }_{{\rm{M}}}$ for large temperatures).

Non-zero mechanical damping reduces the mechanical resonance frequency to ${\omega }_{{\rm{M}}}\sigma$, where $\sigma ={(1-{{\rm{\Gamma }}}^{2}/4{\omega }_{{\rm{M}}}^{2})}^{1/2}\leqslant 1$. It will be convenient to introduce the parameter $\epsilon ={\rm{\Gamma }}/2\sigma {\omega }_{{\rm{M}}}$, which is small in the limit of high mechanical quality (${Q}_{{\rm{M}}}\approx {\omega }_{{\rm{M}}}/{\rm{\Gamma }}\gg 1$). We may solve (8) over a quarter cycle of evolution (appendix B.1), giving

$$\begin{eqnarray*}{X}_{{\rm{M}}} & \to & \sqrt{{\eta }_{{\rm{M}}}}\left[+\epsilon {X}_{{\rm{M}}}+\displaystyle \frac{1}{\sigma }{P}_{{\rm{M}}}+{\rm{\Delta }}{X}_{{\rm{M}}}\right]\\ {P}_{{\rm{M}}} & \to & \sqrt{{\eta }_{{\rm{M}}}}\left[-\epsilon {P}_{{\rm{M}}}-\displaystyle \frac{1}{\sigma }{X}_{{\rm{M}}}+{\rm{\Delta }}{P}_{{\rm{M}}}\right],\end{eqnarray*}$$

where we have introduced noise increments ${\rm{\Delta }}{X}_{{\rm{M}}}$ and ${\rm{\Delta }}{P}_{{\rm{M}}}$ which depend on time integrals of $\hat{\xi }$ (given explicitly in appendix B.1). For a high-quality oscillator ${\eta }_{{\rm{M}}}$ is the amount by which an initial coherent excitation is attenuated over a time of one quarter period (so for an initial excitation amplitude β, ${| \beta | }^{2}\to {\eta }_{{\rm{M}}}{| \beta | }^{2}$ over one quarter cycle); hence, by analogy with ${\eta }_{{\rm{L}}}$, we call this quantity the mechanical efficiency. It is given by ${\eta }_{{\rm{M}}}=\mathrm{exp}\{-\pi \epsilon \}$.

We close this section with some remarks about achieving $\chi =-1$ in the face of experimental imperfections and decoherence. Firstly, we wish to emphasise that many recent pulsed optomechanics experiments have entered a closely related regime where the optomechanical cooperativity (see (13)) exceeds ${\bar{n}}_{{\rm{B}}}$ (e.g. [52, 53]). The essential difference is that our scheme simultaneously requires $\kappa \sim {g}_{0}\sqrt{N}$ to be satisfied, implying that ${g}_{0}\sqrt{N}\gg {\omega }_{{\rm{M}}}$. This necessitates the use of a very large number of pulse photons, because we are by assumption operating with ${g}_{0}\ll {\omega }_{{\rm{M}}}$ and ${\omega }_{{\rm{M}}}\ll \kappa$. Such high-N pulses could be readily generated by acousto-optic or electro-optic (e.g. [36]) modulation of a pump beam over the timescales relevant to our scheme. For example, a pulse lasting 1% of a mechanical period at ${\omega }_{{\rm{M}}}/2\pi =100.2\;{\rm{kHz}}$ requires an average power of ∼10 mW to reach $N\sim {10}^{9}$ (see example experimental parameters given in table 1). The chief technical issue expected to be encountered in this large N limit is heating due to optical absorption, which is often a problem in continuous-wave optomechanics experiments. We present example calculations in appendix B.2 that show, using real experimental parameters, that such heating remains manageable for our pulsed scheme even in the regime of $\chi =-1$.

Table 1.  Exemplary parameters for use with the three-pulse protocol, taken from existing experiments. The mechanical frequency ${\omega }_{{\rm{M}}}$, zero-point coupling rate g₀ and optical linewidth κ are taken from silicon nitride microstring resonators evanescently coupled to silica microsphere whispering gallery modes [54]; the decay rate Γ of comparable oscillators has previously been observed to be as low $2\pi \times 31$ mHz [11, 35]. With these specifications and a bath temperature of 4 K the oscillator may be ground state cooled ($\langle {\hat{n}}_{{\rm{M}}}^{\prime }\rangle =0.606$). Reducing the bath temperature to 50 mK (e.g. [55]) would permit $\langle {\hat{n}}_{{\rm{M}}}^{\prime }\rangle =0.008$ and preparation of the $| 1\rangle$ motional state with infidelity $1-{ \mathcal F }=2.3\%$.

${\omega }_{{\rm{M}}}/2\pi$	100.2 kHz	${\rm{\Gamma }}/2\pi$	31 mHz
$\kappa /2\pi$	25.6 MHz	${g}_{0}/2\pi$	75 Hz
N	$7.28\times {10}^{9}$	C	$7.72\times {10}^{5}$
T	$\left\{\begin{array}{c}4\;{\rm{K}}\\ 5\;\mathrm{mK}\end{array}\right.$	${\bar{n}}_{{\rm{B}}}$	$\left\{\begin{array}{c}8.32\times {10}^{5}\\ 1.04\times {10}^{4}\end{array}\right.$

The three-pulse protocol introduced above may easily be generalised to the case in which the pulses have unequal interaction strengths, which we shall denote by ${\chi }^{(j)}$, where $j=1,2,3$ labels each pulse. This could be arranged by using the control (displacement) pulses to alter the mean number of envelope photons ${N}^{(j)}$ between QND interactions.

It will prove convenient to introduce a set of 'pulse strengths' ${\mu }^{(j)}$ which take into account the level of decoherence present; these are given by

$$\begin{eqnarray*}{\mu }^{(1)} & = & -{\chi }^{(1)}{\mu }^{(2)},\\ {\mu }^{(2)} & = & -\sqrt{{\eta }_{{\rm{L}}}{\eta }_{{\rm{M}}}}{\chi }^{(2)}/\sigma ,\\ {\mu }^{(3)} & = & -{\chi }^{(3)}{\mu }^{(2)}.\end{eqnarray*}$$

Note that in the absence of decoherence the perfect state swap is achieved when ${\mu }^{(j)}=1\ \forall \ j$.

The general transformation between initial and final quadratures is a nonhomogeneous linear map

$$\begin{eqnarray}&&{{\boldsymbol{X}}}^{\prime }=M{\boldsymbol{X}}+{\boldsymbol{F}}\end{eqnarray} \tag{ 9 }$$

in which the vector of quadrature operators ${\boldsymbol{X}}={({X}_{{\rm{M}}},{P}_{{\rm{M}}},{X}_{{\rm{L}}},{P}_{{\rm{L}}})}^{{\rm{T}}}$ is acted upon by the matrix M and there is a vector of zero-mean noise contributions ${\boldsymbol{F}}$ which depends upon ${\mu }^{(2)}$, ${\mu }^{(3)}$ and the noise increments ${\rm{\Delta }}{X}_{{\rm{M}}}$, ${\rm{\Delta }}{P}_{{\rm{M}}}$, X_V and P_V (see appendix C.3).

The explicit matrix form of M is

$$\begin{eqnarray}\left(\begin{array}{cccc}{\mu }^{(1)}-{\eta }_{{\rm{M}}} & 0 & 0 & -{\mu }^{(2)}\\ \epsilon \sigma [{\mu }^{(3)}-{\mu }^{(1)}] & {\mu }^{(3)}-{\eta }_{{\rm{M}}} & {M}_{\mathrm{2,3}} & {\mu }^{(2)}\epsilon \sigma \\ -{\mu }^{(2)}\epsilon \sigma & -{\mu }^{(2)} & {\mu }^{(1)}-{\eta }_{{\rm{L}}} & 0\\ {M}_{\mathrm{3,1}} & 0 & 0 & {\mu }^{(3)}-{\eta }_{{\rm{L}}}\end{array}\right),\end{eqnarray} \tag{ 10 }$$

in which

$$\begin{eqnarray*}{M}_{\mathrm{2,3}} & = & \displaystyle \frac{1}{{\mu }^{(2)}}[{\eta }_{{\rm{L}}}{\mu }^{(3)}+{\mu }^{(1)}({\eta }_{{\rm{M}}}-{\mu }^{(3)})],\\ {M}_{\mathrm{3,1}} & = & \displaystyle \frac{1}{{\mu }^{(2)}}[{\eta }_{{\rm{M}}}{\mu }^{(3)}+{\mu }^{(1)}({\eta }_{{\rm{L}}}-{\mu }^{(3)})].\end{eqnarray*}$$

To achieve a perfect swap the anti-diagonal elements of M are desired to be ±1, while the diagonal blocks—which determine how complete the swap is—should be zero.

In what follows we will restrict ourselves to the case in which the optical and mechanical modes are initially separable and zero-mean ($\langle {\boldsymbol{X}}\rangle ={\bf{0}}$). Note that the latter may always be arranged by local operations on the input states, and that the mean amplitudes are unchanged by the interface.

Ground state cooling is a convenient starting point for several schemes that aim to generate non-classical states of motion [12, 56], fundamental tests of quantum mechanics and many quantum information processing protocols (e.g. see [57]). We shall refer to systems with average occupancies less than one ($\langle {\hat{n}}_{{\rm{M}}}^{\prime }\rangle \lt 1$) as 'ground state cooled' or 'near ground state cooled'. Examples include [13–15].

Clearly the three-pulse state swap protocol introduced here can achieve optomechanical cooling in the bad-cavity limit; a coherent optical pulse carries only vacuum (zero temperature) noise, which can be swapped onto the mechanical system. In order to formulate a ground state cooling criterion we calculate the evolution of the mechanics-light covariance matrix (see appendix D).

The mechanical occupancy after a three-pulse state swap $\langle {\hat{n}}_{{\rm{M}}}^{\prime }\rangle$ (given explicitly in (D3)) may be minimised with respect to the three pulse strengths following the procedure outlined in appendix D. In the high-Q_M limit the minimum is achieved at ${\mu }^{(j)}\approx 1$, yielding

$$\begin{eqnarray}&&{\langle {\hat{n}}_{{\rm{M}}}^{\prime }\rangle }_{\mathrm{min}}\approx \displaystyle \frac{\pi \epsilon }{4}\{3(2{\bar{n}}_{{\rm{B}}}+1)-2{\eta }_{{\rm{L}}}\}+\displaystyle \frac{1}{4}\left\{\displaystyle \frac{1}{{\eta }_{{\rm{L}}}}-1\right\}.\end{eqnarray} \tag{ 11 }$$

If the decoherence is dominated by thermomechanical noise, which is the case in the regime of interest (assuming ${\bar{n}}_{{\rm{B}}}\gg 1$), we may insert (11) into the ground state criterion $\langle {\hat{n}}_{{\rm{M}}}^{\prime }\rangle \lt 1$ to show that ground state cooling is possible only when the additional criterion

$$\begin{eqnarray}&&{\bar{n}}_{{\rm{B}}}\lt \displaystyle \frac{4{\omega }_{{\rm{M}}}}{3\pi {\rm{\Gamma }}}\end{eqnarray} \tag{ 12 }$$

is satisfied. As discussed in section 4 (see table 1), meeting the requirement ${\mu }^{(j)}\approx 1$ is experimentally tractable. The further condition (12) is essentially a requirement of quantum coherent oscillation (i.e. less than one phonon entering the oscillator per mechanical period). Many experimental systems boast Q_M large enough to satisfy this at temperatures achievable with conventional refrigeration technologies (e.g. [11]). Recent devices may permit this regime to even be reached at room temperature (${Q}_{{\rm{M}}}\sim {10}^{8}$ at ${\omega }_{{\rm{M}}}/2\pi \sim 150\;\mathrm{kHz}$) [58].

Other strategies for optomechanical cooling, such as resolved-sideband cooling or cold damping, typically express their criteria for ground state cooling in terms of the optomechanical cooperativity C. We will therefore recast (12) into a requirement on the cooperativity.

C is a dimensionless parameter describes how strongly the quantum back-action of the optomechanical coupling (due to amplitude noise) perturbs the dynamics of the mechanical oscillator: in the steady state

$$\begin{eqnarray*}&&C=\displaystyle \frac{4{g}_{0}^{2}{| \alpha | }^{2}}{{\rm{\Gamma }}\kappa }\ (\mathrm{steady}\;\mathrm{state}).\end{eqnarray*}$$

A natural quantity to use in the pulsed setting is the instantaneous cooperativity time-averaged over the entire three-pulse sequence; thus

$$\begin{eqnarray}C & = & \displaystyle \frac{4{g}_{0}^{2}}{{\rm{\Gamma }}\kappa }\times \displaystyle \frac{\sigma {\omega }_{{\rm{M}}}}{\pi }{\int }_{0}^{\pi /\sigma {\omega }_{{\rm{M}}}}{\rm{d}}t\ {| \alpha (t)| }^{2}\\ & = & \displaystyle \frac{1}{8\pi \epsilon }\left(\displaystyle \frac{{{\mu }^{(1)}}^{2}+{{\mu }^{(3)}}^{2}}{{{\mu }^{(2)}}^{2}}+\displaystyle \frac{{\sigma }^{2}{{\mu }^{(2)}}^{2}}{{\eta }_{{\rm{L}}}{\eta }_{{\rm{M}}}}\right),\end{eqnarray} \tag{ 13 }$$

where in the first line $\alpha (t)$ includes contributions from all three pulses.

Substituting this definition into (12) and noting that ${\mu }^{(j)}\approx 1\ \forall \ j$ for optimum cooling yields

$$\begin{eqnarray}&&{\bar{n}}_{{\rm{B}}}\lt \displaystyle \frac{16}{3(2+{\eta }_{{\rm{L}}}^{-1})}C\leqslant {\left(\displaystyle \frac{4}{3}\right)}^{2}C.\end{eqnarray} \tag{ 14 }$$

Thus, although operating in a pulsed manner and in a fundamentally different regime, our interface's ground state cooling criterion is remarkably similar to that of resolved-sideband cooling (${\bar{n}}_{{\rm{B}}}\lt C$), the current 'gold standard' for optomechanical cooling [42]. We note also that it also resembles the criterion for feedback cooling to near the ground state in the bad-cavity regime, namely ${\bar{n}}_{{\rm{B}}}\lt 8C$ [27]. However, it must be remembered that in (14) the cooperativity must be evaluated with ${\mu }^{(j)}\approx 1$ because a near-perfect state swap is only possible in this regime: merely having a large cooperativity does not guarantee near ground state cooling using the three-pulse protocol.

In the case that the pulse strengths are not equal to unity we find that significant cooling is still possible so long as the deviation in μ is less than a threshold value which scales inversely with $\sqrt{{Q}_{{\rm{M}}}}$ (appendix D.1, figure 5); high-Q_M oscillators are therefore more sensitive to imperfections in μ.

The validity of (14) is confirmed by numerical calculations of the final phonon occupancy for a variety of ${\mu }^{(j)}$, Γ and ${\bar{n}}_{{\rm{B}}}$, as shown in figure 2 (and figure 5). In all figures we have set ${\eta }_{{\rm{L}}}=1$.

Zoom In Zoom Out Reset image size

Figure 2 shows the expected trend of decreasing final occupancy as the oscillator Q_M is increased and the initial temperature is decreased. The analytical threshold for near-ground-state cooling (14) is in excellent agreement with the calculations. Note that the lines of equal heating rate are curved because there is a non-zero incoming noise contribution from the mechanical bath even when the ambient temperature is zero.

Near ground state cooling can be realised using the three-pulse protocol with achievable experimental parameters, as shown in table 1. This cooling will persist for a time proportional to ${({\rm{\Gamma }}{\bar{n}}_{{\rm{B}}})}^{-1}$ after the conclusion of the protocol. We note that these same parameters also permit non-classical state transfer, as discussed below (section 7).

For more complicated quantum states the covariance matrix ceases to be a full description, so we now move to considering the evolution of the system's Wigner function. The Wigner function is a quasiprobability distribution (appendix E) which is permitted to be negative over parts of phase space; negativity is considered to be a 'smoking gun' for non-classical states [56].

Our interface uses only linear interactions and therefore cannot generate Wigner negativity. It can, however, transfer negativity from one system to another. We demonstrate this numerically by considering the preparation of a single phonon Fock state (appendix E.1.2). Fock states have definite energy but completely undefined phase, and as such are highly nonclassical. As may be seen in figure 3 (panels (i)–(iii)), transfer of a single-photon state onto the mechanical oscillator is possible with realistic mechanical parameters (see table 1). Note that the state becomes rotationally asymmetric when the noise entering from the thermal bath is significant, and that increasing the heating rate leads to degradation of the Wigner negativity. This is consistent with thermalisation of the system during the protocol.

Zoom In Zoom Out Reset image size

Higher-number mechanical Fock states may also be prepared. We show the infidelity ($1-{ \mathcal F }$, where ${ \mathcal F }$ is the Schumacher fidelity [59]) of this operation in figure 4 (panel (i)). Note that higher-n Fock states are more difficult to prepare because of their intrinsically larger susceptibility to environmental noise.

Zoom In Zoom Out Reset image size

The interface may be made to perform additional non-trivial operations on the states involved. Here we consider fixing ${\mu }^{(\mathrm{1,3})}=1$ and allowing ${\mu }^{(2)}$ to take on non-unit values. In this case (10) becomes (setting ${\eta }_{{\rm{L}}}={\eta }_{{\rm{M}}}=1$)

$$\begin{eqnarray}M({\mu }^{(2)})=\left(\begin{array}{cccc}0 & 0 & 0 & -{\mu }^{(2)}\\ 0 & 0 & 1/{\mu }^{(2)} & 0\\ 0 & -{\mu }^{(2)} & 0 & 0\\ 1/{\mu }^{(2)} & 0 & 0 & 0\end{array}\right),\end{eqnarray} \tag{ 15 }$$

which may be decomposed into a state transfer, à la (7), followed by local squeezing operations (one on each output) with equal squeezing parameters $\xi =-\mathrm{ln}{\mu }^{(2)}$. Thus ${\mu }^{(2)}\lt 1$ yields position-squeezed output states, and momentum-squeezed outputs may be prepared with ${\mu }^{(2)}\gt 1$.

Zoom In Zoom Out Reset image size

In principle, it is then straightforward to generate unconditionally squeezed motion by transferring the noise of a coherent pulse onto the oscillator with ${\mu }^{(2)}\ne 1$.

A more intriguing application of this squeezed transfer is generation of small-amplitude mechanical Schrödinger cat states (coherent state superpositions), termed 'kitten' states. As has been demonstrated in optical systems, these may be 'bred' into large-amplitude cat states by a number of schemes [60–62]; this could result in a macroscopic superposition of a massive object. Creating such a state is a research goal of foundational interest [12], and is also sufficient to allow universal quantum computation [63, 64] with phonons.

Any state of the form

$$\begin{eqnarray*}&&| {{\rm{\Psi }}}^{(-)}\rangle =\displaystyle \frac{| \alpha \rangle -| -\alpha \rangle }{\sqrt{2(1-{{\rm{e}}}^{-2{| \alpha | }^{2}})}},\end{eqnarray*}$$

is an odd cat state. Taking α to be small and real, we have

$$\begin{eqnarray*}| {{\rm{\Psi }}}^{(-)}\rangle & = & \displaystyle \frac{| \alpha \rangle -| -\alpha \rangle }{\sqrt{2(1-{{\rm{e}}}^{-2{| \alpha | }^{2}})}}=\sqrt{\displaystyle \frac{2}{{{\rm{e}}}^{+{| \alpha | }^{2}}-{{\rm{e}}}^{-{| \alpha | }^{2}}}}\displaystyle \sum _{{\rm{j}}\in \mathrm{odd}}^{\infty }\displaystyle \frac{{\alpha }^{j}}{\sqrt{j!}}| j\rangle \\ & \approx & \displaystyle \frac{\alpha }{| \alpha | }\left(| 1\rangle +\displaystyle \frac{{\alpha }^{2}}{\sqrt{6}}| 3\rangle +{ \mathcal O }({\alpha }^{4})\right).\end{eqnarray*}$$

Note that this contains only odd Fock state contributions.

A squeezer generates correlated pairs of photons; it follows that weakly squeezing $| 1\rangle$ could give an approximation to an odd kitten state (see [65]). For small squeezing parameters ξ

$$\begin{eqnarray*}S(\xi )| 1\rangle & = & \mathrm{exp}\left\{\displaystyle \frac{1}{2}({\xi }^{*}{\hat{a}}^{2}-\xi {{\hat{a}}^{\dagger }}^{2})\right\}| 1\rangle =\displaystyle \sum _{{\rm{j}}=0}^{\infty }{\left(\displaystyle \frac{1}{2}\right)}^{j}\displaystyle \frac{{({\xi }^{*}{\hat{a}}^{2}-\xi {{\hat{a}}^{\dagger }}^{2})}^{j}}{j!}| 1\rangle \\ & \to & | 1\rangle -\xi \sqrt{\displaystyle \frac{3}{2}}| 3\rangle +{ \mathcal O }({\xi }^{2}).\end{eqnarray*}$$

By comparing these we see that the approximate kitten state has an amplitude given by $\xi =-\tfrac{{\alpha }^{2}}{3}$. Thus we obtain ${\mu }^{(2)}=\mathrm{exp}\{{\alpha }^{2}/3\}$.

The Wigner function of such a state is given in appendix E.1.3 and shown with mild decoherence in figure 3 (panel (iv). It is important to note that the Wigner function contains two positive lobes along the ${X}_{{\rm{M}}}^{\prime }$ axis ('alive' and 'dead') separated by a region of quantum interference exhibiting Wigner negativity, proving that coherence between the lobes persists [56]. Again, note that these parameters are experimentally realisable (although requiring very large N to achieve the squeezing shown in figure 3, panel (iv)).

The degree to which the decohering squeezed single phonon state approximates a true cat state is shown in figure 4 (panel (ii)). It is evident that large cat states cannot be produced directly in this manner, as the minimum infidelity rises monotonically with $| \alpha | ;$ this is not due to imperfections in our interface, but from the approximate equality of an odd cat state and a squeezed single phonon state. Environmental noise clearly degrades the preparation process, both through decreasing the purity of the initial mechanical state and through destroying the phonon pair correlations required for squeezing. Even so, it is possible to push the infidelity to low values with realistic parameters (see table 1). We also note that other non-Gaussian quantum states of light can be transferred to the mechanics by our interface thus providing several routes for motional state engineering.

We have shown how open-loop rotations and displacements plus repeated QND interactions may be used to construct an optomechanical interface which operates in the high-bandwidth ($\kappa \gg {\omega }_{{\rm{M}}}$) limit of optomechanics, a natural operating regime for many devices. Our calculations predict that such an interface is within reach of current experimental parameters (see table 1). Furthermore, we have shown how our three-pulse protocol may be extended to performing more complicated operations, including optical–mechanical state transfer with squeezing of the final modes; we note that it is also possible to perform a partial state swap which entangles the final modes, though we have not modelled that process here.

These results show that optomechanical systems in the unresolved sideband limit are capable of performing an optomechanical state-swap, and may prove important in the construction of quantum networks and studies of macroscopic quantum mechanics. Finally, we note that this three-pulse protocol could directly be applied in other circumstances where QND interactions with light are available, such as in spin ensembles.

This work was funded by the Australian Research Council through the Centre of Excellence for Engineered Quantum Systems (EQuS, CE110001013), a Discovery Project grant (MRV, DP140101638) and a Future Fellowship (WPB, FT140100650); LSM further acknowledges support from the Carlsberg Foundation.

The authors thank Nicholas Wyatt, Andreas Næsby, Yauhen 'Eugene' Sachkou and George A. Brawley for useful discussions.

The optomechanical coupling Hamiltonian takes the form given in the main text (see [24, 48]);

$$\begin{eqnarray*}&&\hat{H}={\hslash }{g}_{0}{X}_{{\rm{M}}}{\hat{a}}^{\dagger }\hat{a}.\end{eqnarray*}$$

If we apply a bright pulse to the cavity (on resonance) then we may linearise around the intracavity field's mean amplitude $\alpha (t)$. Doing so amounts to letting $\hat{a}\to \alpha (t)+\delta \hat{a}(t)$ where the displaced annihilation operator $\delta \hat{a}(t)$ is zero mean. Then

$$\begin{eqnarray*}&&\hat{H}={\hslash }{g}_{0}{X}_{{\rm{M}}}({| \alpha (t)| }^{2}+{\alpha }^{*}\delta \hat{a}(t)+\alpha (t)\delta {\hat{a}}^{\dagger }(t)+\delta {\hat{a}}^{\dagger }(t)\delta \hat{a}(t)).\end{eqnarray*}$$

As is common, we neglect the 'small' term $\delta {\hat{a}}^{\dagger }(t)\delta \hat{a}(t)$ because the optomechanical coupling is taken to be weak. Thus we arrive at

$$\begin{eqnarray*}&&\hat{H}={\hslash }| \alpha (t)| {g}_{0}{X}_{{\rm{M}}}{\bar{X}}_{{\rm{L}}}(t)+{\hslash }{g}_{0}{X}_{{\rm{M}}}{| \alpha (t)| }^{2},\end{eqnarray*}$$

as given in the main text. From this we see that the amplitude fluctuations are defined by

$$\begin{eqnarray*}&&{\bar{X}}_{{\rm{L}}}(t)=\displaystyle \frac{\alpha (t)\delta {\hat{a}}^{\dagger }(t)+{\alpha }^{*}(t)\hat{a}(t)}{| \alpha (t)| }.\end{eqnarray*}$$

This reduces to the expression given in the main text if α is real and positive.

The conjugate fluctuations, representing phase noise of the field, may be written

$$\begin{eqnarray*}&&{\bar{P}}_{{\rm{L}}}(t)={\rm{i}}\displaystyle \frac{\alpha (t)\delta {\hat{a}}^{\dagger }(t)-{\alpha }^{*}(t)\hat{a}(t)}{| \alpha (t)| }.\end{eqnarray*}$$

These definitions ensure that the amplitude and phase fluctuations transform correctly when the optical state is displaced. Specifically, if the state is displaced such that $\alpha \to {\rm{i}}\alpha$ then the optical noise is rotated by $\pi /2$. The utility of this definition will become apparent in appendix B.3.

B.1. Free mechanical evolution (rotation)

For a linearly damped oscillator

$$\begin{eqnarray*}{\dot{X}}_{{\rm{M}}} & = & {\omega }_{{\rm{M}}}{P}_{{\rm{M}}},\\ {\dot{P}}_{{\rm{M}}} & = & -{\omega }_{{\rm{M}}}{X}_{{\rm{M}}}+\sqrt{2{\rm{\Gamma }}}\hat{\xi }-{\rm{\Gamma }}{P}_{{\rm{M}}}.\end{eqnarray*}$$

The noise operator $\hat{\xi }$ has the properties [51, 70]

$$\begin{eqnarray*}\langle \hat{\xi }(t)\hat{\xi }({t}^{\prime })\rangle & = & {\int }_{0}^{\infty }{\rm{d}}\omega \ \left\{\displaystyle \frac{\omega }{\pi {\omega }_{{\rm{M}}}}\mathrm{coth}\left(\displaystyle \frac{{\hslash }\omega }{2{k}_{{\rm{B}}}{T}_{{\rm{B}}}}\right)\mathrm{cos}\omega (t-{t}^{\prime })\right\}\\ & & +\displaystyle \frac{{\rm{i}}}{2{\omega }_{{\rm{M}}}}\left(\displaystyle \frac{\partial }{\partial t}-\displaystyle \frac{\partial }{\partial {t}^{\prime }}\right)\delta (t-{t}^{\prime }),\\ [\xi (t),{X}_{{\rm{M}}}(t)] & = & {\rm{i}}\sqrt{2{\rm{\Gamma }}},\\ [Y(t),\xi ({t}^{\prime })] & = & \displaystyle \frac{1}{{\omega }_{{\rm{M}}}}\sqrt{\displaystyle \frac{{\rm{\Gamma }}}{2}}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}^{\prime }}\{u(t-{t}^{\prime })[Y(t),{X}_{{\rm{M}}}({t}^{\prime })]\},\end{eqnarray*}$$

where u is the piecewise step function equal to 1/2 at the origin, 1 for positive arguments and 0 otherwise. $Y(t)$ is an arbitrary system operator. These properties are necessary for preservation of the canonical commutator of X_M and P_M [51]. In the high-temperature limit (where ${\hslash }\omega \ll 2{k}_{{\rm{B}}}{T}_{{\rm{B}}}$ for all relevant frequencies ω) we may approximate the former by

$$\begin{eqnarray*}&&\langle \hat{\xi }(t)\hat{\xi }({t}^{\prime })\rangle =(2{\bar{n}}_{{\rm{B}}}+1)\delta (t-{t}^{\prime })+\displaystyle \frac{{\rm{i}}}{2{\omega }_{{\rm{M}}}}\left(\displaystyle \frac{\partial }{\partial t}-\displaystyle \frac{\partial }{\partial {t}^{\prime }}\right)\delta (t-{t}^{\prime }),\end{eqnarray*}$$

where ${\bar{n}}_{{\rm{B}}}$ is the average occupancy of the mechanical mode when in thermal equilibrium with its bath.

The equations of motion may then be integrated to yield

$$\begin{eqnarray*}{X}_{{\rm{M}}}(t) & = & {{\rm{e}}}^{-{\rm{\Gamma }}t/2}[\mathrm{cos}({\omega }_{{\rm{M}}}\sigma t)+\epsilon \mathrm{sin}({\omega }_{{\rm{M}}}\sigma t)]{X}_{{\rm{M}}}(0)\\ & & +{{\rm{e}}}^{-{\rm{\Gamma }}t/2}\left[\displaystyle \frac{1}{\sigma }\mathrm{sin}({\omega }_{{\rm{M}}}\sigma t){P}_{{\rm{M}}}(0)+{\rm{\Delta }}{X}_{{\rm{M}}}(t)\right],\\ {P}_{{\rm{M}}}(t) & = & {{\rm{e}}}^{-{\rm{\Gamma }}t/2}[\mathrm{cos}({\omega }_{{\rm{M}}}\sigma t)-\epsilon \mathrm{sin}({\omega }_{{\rm{M}}}\sigma t)]{P}_{{\rm{M}}}(0)\\ & & +{{\rm{e}}}^{-{\rm{\Gamma }}t/2}\left[\displaystyle \frac{-1}{\sigma }\mathrm{sin}({\omega }_{{\rm{M}}}\sigma t){X}_{{\rm{M}}}(0)+{\rm{\Delta }}{P}_{{\rm{M}}}(t)\right],\end{eqnarray*}$$

where the damping rate Γ induces an offset to the oscillation frequency, characterised by

$$\begin{eqnarray*}&&\sigma =\sqrt{1-{({\rm{\Gamma }}/2{\omega }_{{\rm{M}}})}^{2}}={(1+{\epsilon }^{2})}^{-1/2},\end{eqnarray*}$$

with $\epsilon ={\rm{\Gamma }}/2{\omega }_{{\rm{M}}}\sigma$ (we assume that the oscillator is underdamped throughout). Noise entering from the thermal bath appears in these expressions as

$$\begin{eqnarray*}{\rm{\Delta }}{X}_{{\rm{M}}}(t) & = & \sqrt{2{\rm{\Gamma }}}{\int }_{0}^{t}{\rm{d}}{t}^{\prime }\ {{\rm{e}}}^{{\rm{\Gamma }}{t}^{\prime }/2}\displaystyle \frac{1}{\sigma }\mathrm{sin}({\omega }_{{\rm{M}}}\sigma (t-{t}^{\prime }))\hat{\xi }({t}^{\prime }),\\ {\rm{\Delta }}{P}_{{\rm{M}}}(t) & = & \sqrt{2{\rm{\Gamma }}}{\int }_{0}^{t}{\rm{d}}{t}^{\prime }\ {{\rm{e}}}^{{\rm{\Gamma }}{t}^{\prime }/2}\hat{\xi }({t}^{\prime })\\ & & \times [\mathrm{cos}({\omega }_{{\rm{M}}}\sigma (t-{t}^{\prime }))-\epsilon \mathrm{sin}({\omega }_{{\rm{M}}}\sigma (t-{t}^{\prime }))].\end{eqnarray*}$$

When evaluated after a quarter of a mechanical cycle (X_M and P_M in the following expressions are implicitly evaluated at time zero unless otherwise stated, while the noise operators are evaluated at the end of the quarter cycle)

$$\begin{eqnarray*}{X}_{{\rm{M}}}(t) & = & \sqrt{{\eta }_{{\rm{M}}}}\left[\epsilon {X}_{{\rm{M}}}+\displaystyle \frac{1}{\sigma }{P}_{{\rm{M}}}+{\rm{\Delta }}{X}_{{\rm{M}}}\right],\\ {P}_{{\rm{M}}}(t) & = & \sqrt{{\eta }_{{\rm{M}}}}\left[-\epsilon {P}_{{\rm{M}}}-\displaystyle \frac{1}{\sigma }{X}_{{\rm{M}}}+{\rm{\Delta }}{P}_{{\rm{M}}}\right]\end{eqnarray*}$$

and the following commutators apply;

$$\begin{eqnarray*}[{\rm{\Delta }}{X}_{{\rm{M}}},{\rm{\Delta }}{P}_{{\rm{M}}}] & = & 2{\rm{i}}({\eta }_{{\rm{M}}}^{-1}-1),\\ [{X}_{{\rm{M}}},{\rm{\Delta }}{X}_{{\rm{M}}}] & = & [{X}_{{\rm{M}}},{\rm{\Delta }}{P}_{{\rm{M}}}]=0,\\ [{P}_{{\rm{M}}},{\rm{\Delta }}{X}_{{\rm{M}}}] & = & +2{\rm{i}}\epsilon ,\\ [{P}_{{\rm{M}}},{\rm{\Delta }}{P}_{{\rm{M}}}] & = & -2{\rm{i}}\sigma {\epsilon }^{2}.\end{eqnarray*}$$

The efficiency ${\eta }_{{\rm{M}}}=\mathrm{exp}\{-\pi \epsilon \}$ is approximately the amount by which a coherent excitation of the oscillator decays over the course of one quarter mechanical period: for an initial coherent excitation amplitude β we have ${| \beta | }^{2}\to {\eta }_{{\rm{M}}}(1/{\sigma }^{2}+{\epsilon }^{2}){| \beta | }^{2}={\eta }_{{\rm{M}}}(1+2{\epsilon }^{2}){| \beta | }^{2}$, which is approximately ${\eta }_{{\rm{M}}}{| \beta | }^{2}$ for high-Q_M oscillators.

B.2. Pulse absorption

Here we show that our scheme does not suffer from significant heating due to absorption (see section 4) by considering an explicit potential realisation. We model a microstring mechanical resonator of thickness ${t}_{{\rm{R}}}=54\;{\rm{nm}}$, width ${w}_{{\rm{R}}}=10$ $\mu {\rm{m}}$ and length ${{\ell }}_{{\rm{R}}}=1\;{\rm{mm}};$ such devices may be fabricated from silicon nitride (see table 1) [54]. Note that all of the parameters below have been realised in existing experiments.

Using the formula for the absorption of a dielectric plate given by [66] (equation (3.7) thereof) and the measured imaginary refractive index of SiN [67] we estimate that the fraction of incident power absorbed by the oscillator is $f\sim {10}^{-5}$. The energy deposited into the resonator by a single pulse is therefore $E\approx {\hslash }{\omega }_{{\rm{p}}}{fN}2{ \mathcal F }/\pi$, where ${ \mathcal F }$ is the cavity finesse (see [68]) and ${\omega }_{{\rm{p}}}$ is the pulse's central frequency. We have taken the energy of a pulse photon to be ${\hslash }{\omega }_{{\rm{p}}}$, ignoring its spectral width.

It takes a characteristic time for thermal energy deposited at the centre of the string to leave, given by ${\omega }_{\mathrm{th}}^{-1}={\rho }_{{\rm{R}}}{c}_{{\rm{R}}}{{\ell }}_{{\rm{R}}}^{2}/{\kappa }_{\mathrm{th}}$, where ${\rho }_{{\rm{R}}}$ is the device's mass density, c_R is its specific heat capacity and ${\kappa }_{\mathrm{th}}$ is its thermal conductivity. We take ${\rho }_{{\rm{R}}}=3.18\times {10}^{3}$ kg m⁻³, ${c}_{{\rm{R}}}=711$ J kg⁻¹ K⁻¹ and ${\kappa }_{\mathrm{th}}=150$ W m⁻¹ K⁻¹ (values for Si₃N₄ taken from [69], tables 40, 88 and 104). Given the geometry described here, ${\omega }_{{\rm{M}}}\gg {\omega }_{\mathrm{th}}\approx 260\;{\rm{Hz}}$. This means that we may treat the dynamics of the heat distribution as very slow compared to the mechanical motion of interest.

During the full protocol time of half a mechanical period, thermal energy deposited at the centre of the string has propagated a distance (root-mean-squared) of $L=\sqrt{2\pi {\kappa }_{\mathrm{th}}/{\rho }_{{\rm{R}}}{c}_{{\rm{R}}}{\omega }_{{\rm{M}}}}$. Thus the volume of heated material may be approximated as ${V}_{{\rm{R}}}\approx {t}_{{\rm{R}}}{w}_{{\rm{R}}}L$. The resulting change in temperature of the material over this volume is $\delta T=E/\pi {\rho }_{{\rm{R}}}{c}_{{\rm{R}}}{V}_{{\rm{R}}}$. As a consequence, the occupancy of the bath seen by the mechanical mode changes by (at worst)

$$\begin{eqnarray*}&&\delta {\bar{n}}_{{\rm{B}}}=3{k}_{{\rm{B}}}\delta T/{\hslash }{\omega }_{{\rm{M}}}=3\displaystyle \frac{2{{fk}}_{{\rm{B}}}{\omega }_{{\rm{p}}}{ \mathcal F }}{{\pi }^{2}{\omega }_{{\rm{M}}}{\rho }_{{\rm{R}}}{c}_{{\rm{R}}}{V}_{{\rm{R}}}}N.\end{eqnarray*}$$

The overall factor of three accounts for the fact that there are three pulses.

If we employ a pulse with photons at a central wavelength of 1559 nm [54] we obtain $\delta {\bar{n}}_{{\rm{B}}}\approx {10}^{-9}\times { \mathcal F }N$. Thus for a finesse of 100 and photon number as in table 1 the change in the effective bath occupancy is approximately 10% at 50 mK.

This result is in many ways a worst-case scenario. Firstly, the change in temperature seen by the mode of interest is likely to be lower than calculated here because the mode extends over a larger region of the string. Secondly, alternative resonator geometries (which are not long and thin) can permit much lower $\delta {\bar{n}}_{{\rm{B}}}$ because of their larger volume and faster thermal diffusion times. We are therefore confident that heating will not (in principle) prevent one from achieving $\chi =-1$.

B.3. Optical displacements and loss

A displacement of the optical field in any direction not parallel to its initial phase angle has the effect of rotating the quadrature fluctuations (because the instantaneous fluctuations have been defined to be parallel and perpendicular to the mean amplitude, see appendix A). For instance, consider the special case of a displacement sending the mean amplitude α to ${\rm{i}}\alpha ;$ this yields the transformation

$$\begin{eqnarray*}{X}_{{\rm{L}}} & \to & {P}_{{\rm{L}}},\\ {P}_{{\rm{L}}} & \to & -{X}_{{\rm{L}}},\end{eqnarray*}$$

which is precisely a $\pi /2$ rotation of the noise (by contrast, a rotation of the optical field leaves the noise unchanged, because the fluctuations are also rotated).

Such an operation may be realised by interfering the circulating pulse with a very bright incoming pulse, temporally matched to the former but offset in phase, on a highly asymmetric beamsplitter. Given that the noise associated with this operation may be made arbitrarily small in principle (by selecting an extreme asymmetry and applying a more powerful displacement pulse) we will neglect it, or else absorb it into the total efficiency ${\eta }_{{\rm{L}}}$ defined below.

Any optical losses due to inefficiencies within the recirculating loop may be partially compensated by altering the magnitude and phase angle of the displacement pulse, but contamination of the fluctuation quadratures by vacuum noise cannot be reversed. Losses may be modelled as an asymmetric beamsplitter which transmits $\sqrt{{\eta }_{{\rm{L}}}}$ of the circulating pulse ($0\lt {\eta }_{{\rm{L}}}\leqslant 1$) and couples in vacuum noise proportional to $\sqrt{1-{\eta }_{{\rm{L}}}}$.

$$\begin{eqnarray*}{X}_{{\rm{L}}} & \to & \sqrt{{\eta }_{{\rm{L}}}}{X}_{{\rm{L}}}+\sqrt{1-{\eta }_{{\rm{L}}}}{X}_{{\rm{V}}},\\ {P}_{{\rm{L}}} & \to & \sqrt{{\eta }_{{\rm{L}}}}{P}_{{\rm{L}}}+\sqrt{1-{\eta }_{{\rm{L}}}}{P}_{{\rm{V}}}.\end{eqnarray*}$$

C.1. QND interaction

The QND interaction (including open loop displacement of the mechanics) may be written as

$$\begin{eqnarray*}&&{\boldsymbol{X}}\to {M}_{\mathrm{QND}}^{(j)}{\boldsymbol{X}}\end{eqnarray*}$$

with

$$\begin{eqnarray*}{M}_{\mathrm{QND}}^{(j)}=\left(\begin{array}{cccc}1 & 0 & 0 & 0\\ 0 & 1 & {\chi }^{(j)} & 0\\ 0 & 0 & 1 & 0\\ {\chi }^{(j)} & 0 & 0 & 1\end{array}\right).\end{eqnarray*}$$

C.2. Evolution during pulse-delay-rotate cycle

We now consider the net evolution induced by a QND interaction followed by a quarter-cycle delay and an optical displacement. The full protocol consists of two of these cycles followed by a further QND interaction.

The evolution induced by a single pulse-displace-delay cycle may be expressed as a nonhomogeneous linear equation

$$\begin{eqnarray*}&&{\boldsymbol{X}}\to {M}^{(j)}{\boldsymbol{X}}+{{\boldsymbol{F}}}^{(j)},\end{eqnarray*}$$

where ${M}^{(j)}$ is a square matrix characterised by the interaction strength ${\chi }^{(j)}=\tfrac{-8{g}_{0}\sqrt{{N}^{(j)}}}{\kappa }$ and the efficiencies ${\eta }_{{\rm{L}}}$ and ${\eta }_{{\rm{M}}};$

$$\begin{eqnarray*}{M}^{(j)}=\left(\begin{array}{cccc}\epsilon \sqrt{{\eta }_{{\rm{M}}}} & \sqrt{{\eta }_{{\rm{M}}}}/\sigma & \sqrt{{\eta }_{{\rm{M}}}}{\chi }^{(j)}/\sigma & 0\\ -\sqrt{{\eta }_{{\rm{M}}}}/\sigma & -\epsilon \sqrt{{\eta }_{{\rm{M}}}} & -\epsilon {\chi }^{(j)}\sqrt{{\eta }_{{\rm{M}}}} & 0\\ \sqrt{{\eta }_{{\rm{L}}}}{\chi }^{(j)} & 0 & 0 & \sqrt{{\eta }_{{\rm{L}}}}\\ 0 & 0 & -\sqrt{{\eta }_{{\rm{L}}}} & 0\end{array}\right).\end{eqnarray*}$$

The inhomogeneous term corresponds to thermal and vacuum noises; ${{\boldsymbol{F}}}^{(j)}={{\boldsymbol{f}}}^{(j)}+(1-{\eta }_{{\rm{L}}}){{\boldsymbol{v}}}^{(j)}$, where

$$\begin{eqnarray*}&&{{\boldsymbol{f}}}^{(j)}=\sqrt{{\eta }_{{\rm{M}}}}{({\rm{\Delta }}{X}_{{\rm{M}}}^{(j)},{\rm{\Delta }}{P}_{{\rm{M}}}^{(j)},0,0)}^{{\rm{T}}}\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&{{\boldsymbol{v}}}^{(j)}={(0,0,{X}_{{\rm{V}}}^{(j)},{P}_{{\rm{V}}}^{(j)})}^{{\rm{T}}}\end{eqnarray*}$$

in which ${X}_{{\rm{V}}}^{(j)}$ (${P}_{{\rm{V}}}^{(j)}$) are vacuum noise operators and ${\rm{\Delta }}{P}_{{\rm{M}}}^{(j)}$ and ${\rm{\Delta }}{X}_{{\rm{M}}}^{(j)}$ are evaluated at $t=\pi /2{\omega }_{{\rm{M}}}\sigma$ (i.e. at one-quarter mechanical period).

The noise has the following covariance, with other terms vanishing;

$$\begin{eqnarray*}\langle {\rm{\Delta }}{X}_{{\rm{M}}}^{(j)}{\rm{\Delta }}{X}_{{\rm{M}}}^{(k)}\rangle & = & \langle {\rm{\Delta }}{P}_{{\rm{M}}}^{(j)}{\rm{\Delta }}{P}_{{\rm{M}}}^{(k)}\rangle ={\delta }_{{\rm{j}},{\rm{k}}}(2{\bar{n}}_{{\rm{B}}}+1)\left(\displaystyle \frac{1}{{\eta }_{{\rm{M}}}}-1-2{\epsilon }^{2}\right),\\ {\langle {\rm{\Delta }}{X}_{{\rm{M}}}^{(j)}{\rm{\Delta }}{P}_{{\rm{M}}}^{(k)}\rangle }_{{\rm{S}}} & = & 4{\delta }_{{\rm{j}},{\rm{k}}}(2{\bar{n}}_{{\rm{B}}}+1)\displaystyle \frac{\epsilon }{\sigma },\\ \langle {X}_{{\rm{V}}}^{(j)}{X}_{{\rm{V}}}^{(k)}\rangle & = & \langle {P}_{{\rm{V}}}^{(j)}{P}_{{\rm{V}}}^{(k)}\rangle ={\delta }_{{\rm{j}},{\rm{k}}},\\ {\langle {X}_{{\rm{V}}}^{(j)}{P}_{{\rm{V}}}^{(k)}\rangle }_{{\rm{S}}} & = & 0.\end{eqnarray*}$$

${\delta }_{{\rm{j}},{\rm{k}}}$ is the Kronecker delta.

C.3. Full evolution

Taking the evolution given above and iterating (remembering that the final pulsed interaction is not followed by a rotation and optical displacement) yields

$$\begin{eqnarray*}&&{{\boldsymbol{X}}}^{\prime }=M{\boldsymbol{X}}+{\boldsymbol{F}},\end{eqnarray*}$$

where

$$\begin{eqnarray*}M & = & {M}_{\mathrm{QND}}^{(3)}\;{M}^{(2)}{M}^{(1)},\\ \;{\boldsymbol{F}} & = & {M}_{\mathrm{QND}}^{(3)}[{M}^{(2)}{{\boldsymbol{F}}}^{(1)}+{{\boldsymbol{F}}}^{(2)}].\end{eqnarray*}$$

The explicit matrix form of M is

$$\begin{eqnarray*}\left(\begin{array}{cccc}{\mu }^{(1)}-{\eta }_{{\rm{M}}} & 0 & 0 & -{\mu }^{(2)}\\ \epsilon \sigma [{\mu }^{(3)}-{\mu }^{(1)}] & {\mu }^{(3)}-{\eta }_{{\rm{M}}} & {M}_{\mathrm{2,3}} & {\mu }^{(2)}\epsilon \sigma \\ -{\mu }^{(2)}\epsilon \sigma & -{\mu }^{(2)} & {\mu }^{(1)}-{\eta }_{{\rm{L}}} & 0\\ {M}_{\mathrm{3,1}} & 0 & 0 & {\mu }^{(3)}-{\eta }_{{\rm{L}}}\end{array}\right),\end{eqnarray*}$$

in which

$$\begin{eqnarray*}{M}_{\mathrm{2,3}} & = & \displaystyle \frac{1}{{\mu }^{(2)}}[{\eta }_{{\rm{L}}}{\mu }^{(3)}+{\mu }^{(1)}({\eta }_{{\rm{M}}}-{\mu }^{(3)})],\\ {M}_{\mathrm{3,1}} & = & \displaystyle \frac{1}{{\mu }^{(2)}}[{\eta }_{{\rm{M}}}{\mu }^{(3)}+{\mu }^{(1)}({\eta }_{{\rm{L}}}-{\mu }^{(3)})].\end{eqnarray*}$$

It is instructive to consider the form of M in some ideal cases.

C.3.1. State swap

Setting ${\mu }^{(j)}=1$ with no decoherence ($\epsilon =0$, ${\eta }_{{\rm{M}}}={\eta }_{{\rm{L}}}=1$) yields

$$\begin{eqnarray*}M=\left(\begin{array}{cccc}0 & 0 & 0 & -1\\ 0 & 0 & 1 & 0\\ 0 & -1 & 0 & 0\\ 1 & 0 & 0 & 0\end{array}\right).\end{eqnarray*}$$

C.3.2. Squeezed state swap

Setting ${\mu }^{(\mathrm{1,3})}=1$ with no decoherence gives

$$\begin{eqnarray*}M=\left(\begin{array}{cccc}0 & 0 & 0 & -{\mu }^{(2)}\\ 0 & 0 & 1/{\mu }^{(2)} & 0\\ 0 & -{\mu }^{(2)} & 0 & 0\\ 1/{\mu }^{(2)} & 0 & 0 & 0\end{array}\right).\end{eqnarray*}$$

This may be factorised into a state swap (as given above) followed by the application of the two local squeezing operations

$$\begin{eqnarray*}S & = & \mathrm{exp}\left\{\displaystyle \frac{1}{2}{\xi }^{*}\;{\hat{a}}^{2}-\displaystyle \frac{1}{2}\xi {{\hat{a}}^{\dagger }}^{2}\right\}\mathrm{exp}\left\{\displaystyle \frac{1}{2}{\xi }^{*}\;{\hat{b}}^{2}-\displaystyle \frac{1}{2}\xi \;{{\hat{b}}^{\dagger }}^{2}\right\}\\ & \to & \left(\begin{array}{cccc}{\mu }^{(2)} & 0 & 0 & 0\\ 0 & 1/{\mu }^{(2)} & 0 & 0\\ 0 & 0 & {\mu }^{(2)} & 0\\ 0 & 0 & 0 & 1/{\mu }^{(2)}\end{array}\right),\end{eqnarray*}$$

which have equal squeezing parameters $\xi =-\mathrm{ln}{\mu }^{(2)}$.

Under the assumptions listed in the main text (section 6), plus the assumption that the noise is uncorrelated with the initial state, it is straightforward to show that the output covariance matrix

$$\begin{eqnarray*}&&{V}_{{\boldsymbol{X}}{\boldsymbol{X}}}^{\prime }={\mathfrak{R}}\{\langle {{\boldsymbol{X}}}^{\prime }{{{\boldsymbol{X}}}^{\prime }}^{{\rm{T}}}\rangle \}\end{eqnarray*}$$

is related to the input by (see (9))

$$\begin{eqnarray*}&&{V}_{{\boldsymbol{X}}{\boldsymbol{X}}}^{\prime }={{MV}}_{{\boldsymbol{X}}{\boldsymbol{X}}}{M}^{{\rm{T}}}+{V}_{{\boldsymbol{F}}{\boldsymbol{F}}}.\end{eqnarray*}$$

The first diagonal block of ${V}_{{\boldsymbol{X}}{\boldsymbol{X}}}^{\prime }$ corresponds to the output mechanical covariance matrix, the second represents the output optical state, and the antidiagonal blocks describe correlations between the two subsystems.

In the high-Q_M limit we may expand this result about $\epsilon =0$ and determine the final mechanical variances

$$\begin{eqnarray}\langle {{X}_{{\rm{M}}}^{\prime }}^{2}\rangle & = & \langle {X}_{{\rm{M}}}^{2}\rangle ({\mu }^{(1)}-1)({\mu }^{(1)}+2\pi \epsilon -1)+\langle {P}_{{\rm{L}}}^{2}\rangle {{\mu }^{(2)}}^{2}\\ & & +2\pi \epsilon (2{\bar{n}}_{{\rm{B}}}+1)+\left(\displaystyle \frac{1}{{\eta }_{{\rm{L}}}}-1\right){{\mu }^{(2)}}^{2},\end{eqnarray} \tag{ D1a }$$

$$\begin{eqnarray}\langle {{P}_{{\rm{M}}}^{\prime }}^{2}\rangle & = & \langle {P}_{{\rm{M}}}^{2}\rangle ({\mu }^{(3)}-1)({\mu }^{(3)}+2\pi \epsilon -1)+{\langle {X}_{{\rm{M}}}{P}_{{\rm{M}}}\rangle }_{{\rm{S}}}2\epsilon ({\mu }^{(3)}-1)({\mu }^{(3)}-{\mu }^{(1)})\\ & & +{\langle {X}_{{\rm{L}}}{P}_{{\rm{L}}}\rangle }_{{\rm{S}}}2\epsilon ({\mu }^{(1)}+{\mu }^{(3)}{\eta }_{{\rm{L}}}-{\mu }^{(1)}{\mu }^{(3)})+\langle {X}_{{\rm{L}}}^{2}\rangle {{\mu }^{(2)}}^{-2}\{{\mu }^{(1)}({\mu }^{(3)}-1)-{\eta }_{{\rm{L}}}{\mu }^{(3)}\}\\ & & \times \{{\mu }^{(1)}({\mu }^{(3)}+2\pi \epsilon -1)-{\eta }_{{\rm{L}}}{\mu }^{(3)}\}+\pi \epsilon (2{\bar{n}}_{{\rm{B}}}+1)\{2-{\mu }^{(3)}(2-{\mu }^{(3)})\}\\ & & +{\left(\displaystyle \frac{{\mu }^{(3)}}{{\mu }^{(2)}}\right)}^{2}(1-{\eta }_{{\rm{L}}}^{2}).\end{eqnarray} \tag{ D1b }$$

It is safe to assume that the initial mechanical state is thermal, viz

$$\begin{eqnarray*}&&\langle {X}_{{\rm{M}}}^{2}\rangle =\langle {P}_{{\rm{M}}}^{2}\rangle =2{\bar{n}}_{{\rm{B}}}+1,\ {\langle {X}_{{\rm{M}}}{P}_{{\rm{M}}}\rangle }_{{\rm{S}}}=0\end{eqnarray*}$$

and because a coherent pulse carries only vacuum noise we have

$$\begin{eqnarray*}&&\langle {X}_{{\rm{L}}}^{2}\rangle =\langle {P}_{{\rm{L}}}^{2}\rangle =1,\ {\langle {X}_{{\rm{L}}}{P}_{{\rm{L}}}\rangle }_{{\rm{S}}}=0.\end{eqnarray*}$$

The final occupancy of the mechanical oscillator is then readily calculated using

$$\begin{eqnarray}&&\langle {\hat{n}}_{{\rm{M}}}^{\prime }\rangle =\displaystyle \frac{1}{4}[\langle {{X}_{{\rm{M}}}^{\prime }}^{2}\rangle +\langle {{P}_{{\rm{M}}}^{\prime }}^{2}\rangle -2].\end{eqnarray} \tag{ D2 }$$

This yields

$$\begin{eqnarray}4\langle {\hat{n}}_{{\rm{M}}}^{\prime }\rangle & = & (2{\bar{n}}_{{\rm{B}}}+1)[({\mu }^{(1)}-1)({\mu }^{(1)}+2\pi \epsilon -1)\\ & & +({\mu }^{(3)}-1)({\mu }^{(3)}+2\pi \epsilon -1)\\ & & +\pi \epsilon (4-{\mu }^{(3)}(2-{\mu }^{(3)}))]\\ & & +\displaystyle \frac{1}{{{\mu }^{(2)}}^{2}}[{{\mu }^{(3)}}^{2}-2{\eta }_{{\rm{L}}}{\mu }^{(1)}{\mu }^{(3)}({\mu }^{(3)}+\pi \epsilon -1)\\ & & +{{\mu }^{(1)}}^{2}({\mu }^{(3)}-1)({\mu }^{(3)}+2\pi \epsilon -1)]\\ & & +\displaystyle \frac{{{\mu }^{(2)}}^{2}}{{\eta }_{{\rm{L}}}}-2.\end{eqnarray} \tag{ D3 }$$

To minimise this, note that $\langle {\hat{n}}_{{\rm{M}}}^{\prime }\rangle$ may be divided into a term proportional to $(2{\bar{n}}_{{\rm{B}}}+1)$ and a term independent of the thermal bath. The former is dependent only upon ${\mu }^{(1)}$ and ${\mu }^{(3)}$, while the latter depends on all three strengths. We therefore choose to approximately minimise the occupancy by minimising the two terms separately.

Minimising the bath-dependent term yields, to first order in

$$\begin{eqnarray*}&&{\mu }^{(\mathrm{1,3})}\approx 1-\pi \epsilon \approx {\eta }_{{\rm{M}}}.\end{eqnarray*}$$

We then substitute these into the remaining noise contribution and select ${\mu }^{(2)}$ accordingly. Thus for optimum cooling we require

$$\begin{eqnarray*}{\mu }^{(1)} & \approx & {\eta }_{{\rm{M}}}\\ {\mu }^{(2)} & \approx & \displaystyle \frac{1+{\eta }_{{\rm{M}}}}{2}{\eta }_{{\rm{L}}}^{1/4}\\ {\mu }^{(3)} & \approx & {\eta }_{{\rm{M}}}.\end{eqnarray*}$$

As ${\eta }_{{\rm{M}}}$ and ${\eta }_{{\rm{L}}}$ approach 1 we obtain ${\mu }^{(j)}=1$ to a very good approximation, as stated previously.

Evaluated at ${\mu }^{(j)}=1$

$$\begin{eqnarray*}&&{\langle {\hat{n}}_{{\rm{M}}}^{\prime }\rangle }_{\mathrm{min}}=\displaystyle \frac{\pi \epsilon }{4}\{3(2{\bar{n}}_{{\rm{B}}}+1)-2{\eta }_{{\rm{L}}}\}+\displaystyle \frac{1}{4}\left\{\displaystyle \frac{1}{{\eta }_{{\rm{L}}}}-1\right\}.\end{eqnarray*}$$

This result is given as (11) in the main text.

Note that in order to ground state cool one must satisfy ${\eta }_{{\rm{L}}}\gt 1/5$ in addition to the criteria given in the main text (14).

D.1. Tolerance to pulse strength variations

It is clear from figure 5 that the system becomes less tolerant to variations in ${\mu }^{(j)}$ for very high-Q_M oscillators, which means that the mean number of photons in the coherent envelope of each pulse must be controlled to a very high degree when attempting to reach the mechanical ground state.

To quantify this we consider the case where all pulse strengths are equal (denoted by μ with no superscript). We first find the two values of μ for which

$$\begin{eqnarray*}&&\langle {\hat{n}}_{{\rm{M}}}^{\prime }\rangle =2{\langle {\hat{n}}_{{\rm{M}}}^{\prime }\rangle }_{\mathrm{min}}.\end{eqnarray*}$$

The difference of these values is used as a proxy for the width of the cooling dip (see figure 5, main text). If ${\bar{n}}_{{\rm{B}}}\gg 1$ and ${\eta }_{{\rm{L}}}\gtrsim 0.1$ this reduces to

$$\begin{eqnarray*}&&\mathrm{width}\approx \sqrt{6\pi \epsilon }.\end{eqnarray*}$$

Thus the tolerance of the protocol to imperfections in the pulse strengths scales as ${\epsilon }^{1/2}$. This is a very strong dependence in the high-Q_M limit; for instance, a mechanical Q_M of 10⁶ would require the value of μ to be stabilised to approximately one part in a thousand.

To translate this into terms of photon number, note that width $\approx {\rm{\Delta }}\mu \approx -{\rm{\Delta }}\chi$, where ${\rm{\Delta }}\mu$ (${\rm{\Delta }}\chi$) is the tolerance of the protocol to μ (χ). Then we may expand the definition of χ to see that

$$\begin{eqnarray*}&&{\rm{\Delta }}\chi \approx \displaystyle \frac{-8{g}_{0}}{\kappa }\sqrt{N}\left(\displaystyle \frac{{\rm{\Delta }}N}{2N}\right),\end{eqnarray*}$$

which becomes ${\rm{\Delta }}\chi \approx -{\rm{\Delta }}N/2N$ in the case that $\chi \approx -1$. Thus ${\rm{\Delta }}N/N\approx 2\sqrt{6\pi \epsilon }$ is the tolerable fractional variation of the pulse photon number amplitude. For Q_M of 10⁶ the fractional variation needs to be less than approximately $0.6\%$.

This $\sqrt{\epsilon }$ scaling may be explained by considering the two ways in which the thermal bath contributes to the final mechanical temperature. Firstly, there is noise entering the system during the protocol; this component is reasonably weakly dependent on μ, and tends to set the minimum attainable final temperature. The second contribution is from the initial thermal state of the oscillator. The amount of this noise which is swapped out of the mechanical system depends strongly on μ, and tends to set the shape of the $\langle {\hat{n}}_{{\rm{M}}}^{\prime }(\mu )\rangle$ curve (see figure 5, main text). If the former dominates then the strong μ dependence of the second term is washed out.

When considering non-Gaussian input states, such as Fock states, the covariance matrix is no longer an adequate description of the state. We therefore consider the Wigner function description of quantum states; for n modes this is defined by

$$\begin{eqnarray}&&W({\boldsymbol{r}})=\displaystyle \frac{1}{{(2\pi )}^{2n}}\int {{\rm{d}}}^{2n}{\boldsymbol{\beta }}\ {\rm{Tr}}\{\hat{\rho }\hat{D}({\boldsymbol{\beta }})\}\mathrm{exp}\{-{\rm{i}}{\boldsymbol{r}}\cdot {\rm{\Omega }}{\boldsymbol{\beta }}\},\end{eqnarray} \tag{ E1 }$$

in which the $2n\times 2n$ matrix Ω is the n-fold block diagonal matrix with diagonal blocks

$$\begin{eqnarray*}\varpi =\left(\begin{array}{rc}0 & 1\\ -1 & 0\end{array}\right)\end{eqnarray*}$$

and all other entries equal to zero (note that ${{\rm{\Omega }}}^{{\rm{T}}}=-{\rm{\Omega }}$ and ${\rm{\Omega }}{{\rm{\Omega }}}^{{\rm{T}}}={{\rm{\Omega }}}^{{\rm{T}}}{\rm{\Omega }}=-{\mathfrak{1}}$, with ${\mathfrak{1}}$ being the identity matrix) [50]. Points in phase space are denoted ${{\boldsymbol{r}}}^{{\rm{T}}}=({x}_{1},{p}_{1},{x}_{2},\ldots )$ and we have suppressed the bounds of integration, which are at $\pm \infty$ for each quadrature. $\hat{D}({\boldsymbol{\beta }})$ is the multi-mode displacement operator with ${{\boldsymbol{\beta }}}^{{\rm{T}}}=({\mathfrak{R}}\{{\beta }_{1}\},{\mathfrak{I}}\{{\beta }_{1}\},{\mathfrak{R}}\{{\beta }_{2}\},\ldots )$ giving the displacement of each mode.

The Wigner function is Fourier dual to the characteristic function $\chi ({\boldsymbol{\beta }})={\rm{Tr}}\{\hat{\rho }\hat{D}({\boldsymbol{\beta }})\}$, viz

$$\begin{eqnarray}&&\chi ({\boldsymbol{\beta }})=\int {{\rm{d}}}^{2n}{\boldsymbol{r}}\ W({\boldsymbol{r}})\mathrm{exp}\{+{\rm{i}}{\boldsymbol{r}}\cdot {\rm{\Omega }}{\boldsymbol{\beta }}\}\end{eqnarray} \tag{ E2 }$$

and knowledge of either is sufficient to reconstruct the density operator $\hat{\rho }$ [50]. It will prove useful to work in the reciprocal domain (${\boldsymbol{\beta }}$) for the calculations which follow.

E.1. Wigner functions used in this work

E.1.1. Gaussian states

Gaussian Wigner functions have the following form, where $\delta {\boldsymbol{r}}={\boldsymbol{r}}-\langle {\boldsymbol{X}}\rangle ={({x}_{1}-\langle {X}_{1}\rangle ,{p}_{1}-\langle {P}_{1}\rangle ,\ldots )}^{{\rm{T}}}$:

$$\begin{eqnarray*}W({\boldsymbol{r}}) & = & \displaystyle \frac{1}{{(2\pi )}^{n}\sqrt{| V| }}\mathrm{exp}\left\{-\displaystyle \frac{1}{2}{\delta {\boldsymbol{r}}}^{{\rm{T}}}{V}^{-1}\delta {\boldsymbol{r}}\right\}.\\ \chi ({\boldsymbol{\beta }}) & = & \mathrm{exp}\left\{-\displaystyle \frac{1}{2}{{\boldsymbol{\beta }}}^{{\rm{T}}}{{\rm{\Omega }}}^{{\rm{T}}}V{\rm{\Omega }}{\boldsymbol{\beta }}\right\}{{\rm{e}}}^{{\rm{i}}\langle {\boldsymbol{X}}\rangle \cdot {\rm{\Omega }}{\boldsymbol{\beta }}}.\end{eqnarray*}$$

They are entirely specified by the mean values and $\langle {\boldsymbol{X}}\rangle$ and the covariance matrix $V={\mathfrak{R}}\{\langle {\boldsymbol{X}}{{\boldsymbol{X}}}^{{\rm{T}}}\rangle \}$.

E.1.2. Fock states

The Fock states $| n\rangle$ are non-Gaussian, having Wigner functions given by the Laguerre polynomials ${{ \mathcal L }}_{n}$. We give these for a single mode only

$$\begin{eqnarray*}{\chi }_{n}(\beta ) & = & {{\rm{e}}}^{-\displaystyle \frac{1}{2}{| \beta | }^{2}}{{ \mathcal L }}_{n}({| \beta | }^{2}),\\ {W}_{n}({\boldsymbol{r}}) & = & \displaystyle \frac{1}{2\pi }{(-1)}^{n}\mathrm{exp}\left\{-\displaystyle \frac{1}{2}{\boldsymbol{r}}\cdot {\boldsymbol{r}}\right\}{{ \mathcal L }}_{n}({\boldsymbol{r}}\cdot {\boldsymbol{r}}).\end{eqnarray*}$$

E.1.3. Cat states

The final class of non-Gaussian state treated in this paper are the odd (−) cat states. For completeness we include the expressions for the even (+) cat states too. Given a finite amplitude α (for a single mode)

$$\begin{eqnarray*}{\chi }_{\alpha }^{(\pm )}(\beta ) & = & {(1\pm {{\rm{e}}}^{-2{| \alpha | }^{2}})}^{-1}{{\rm{e}}}^{-\displaystyle \frac{1}{2}{| \beta | }^{2}}[\mathrm{cosh}[{\alpha }^{*}\beta -{\beta }^{*}\alpha ]\pm {{\rm{e}}}^{-2{| \alpha | }^{2}}\mathrm{cosh}[{\alpha }^{*}\beta +{\beta }^{*}\alpha ]],\\ {W}_{\alpha }^{(\pm )}({\boldsymbol{r}}) & = & \displaystyle \frac{{{\rm{e}}}^{-\tfrac{1}{2}{\boldsymbol{r}}\cdot {\boldsymbol{r}}}}{2\pi (1\pm {{\rm{e}}}^{-2{| \alpha | }^{2}})}[{{\rm{e}}}^{-2{| \alpha | }^{2}}\mathrm{cosh}(2{\boldsymbol{r}}\cdot {\boldsymbol{\alpha }})\pm \mathrm{cos}(2{\boldsymbol{r}}\cdot \varpi {\boldsymbol{\alpha }})].\end{eqnarray*}$$

E.2. Lossy evolution

The output characteristic function may be written (according to its definition in the Heisenberg picture, for a linear system) as

$$\begin{eqnarray*}{\chi }^{\prime }({\boldsymbol{\beta }}) & = & {\rm{Tr}}\{\hat{\rho }\mathrm{exp}[{\rm{i}}{{\boldsymbol{X}}}^{\prime }\cdot {\rm{\Omega }}{\boldsymbol{\beta }}]\}\\ & = & {\rm{Tr}}\{\hat{\rho }\mathrm{exp}[{\rm{i}}(M{\boldsymbol{X}}+{\boldsymbol{F}})\cdot {\rm{\Omega }}{\boldsymbol{\beta }}]\}\\ & = & {\rm{Tr}}\{\hat{\rho }\mathrm{exp}[{\rm{i}}{{\boldsymbol{\beta }}}^{{\rm{T}}}{{\rm{\Omega }}}^{{\rm{T}}}M{\boldsymbol{X}}]\mathrm{exp}[{\rm{i}}{\boldsymbol{F}}\cdot {\rm{\Omega }}{\boldsymbol{\beta }}]\}\\ & = & {\rm{Tr}}\{\hat{\rho }\mathrm{exp}[{\rm{i}}{{\boldsymbol{\beta }}}^{{\rm{T}}}{{\rm{\Omega }}}^{{\rm{T}}}M{\boldsymbol{X}}]\}{\rm{Tr}}\{\hat{\rho }\mathrm{exp}[{\rm{i}}{\boldsymbol{F}}\cdot {\rm{\Omega }}{\boldsymbol{\beta }}]\}.\end{eqnarray*}$$

In going from line two to line three we used the Baker–Campbell–Haussdorf lemma, and to arrive at line four we noted that expectation values of system and noise operators factorise under the approximations used here.

If we define a vector ${\boldsymbol{\gamma }}={\rm{\Omega }}{M}^{{\rm{T}}}{\rm{\Omega }}{\boldsymbol{\beta }}$ then it is clear that the exponent of the first trace above becomes

$$\begin{eqnarray*}{{\boldsymbol{\beta }}}^{{\rm{T}}}{{\rm{\Omega }}}^{{\rm{T}}}M{\boldsymbol{X}} & = & -{{\boldsymbol{\beta }}}^{{\rm{T}}}{{\rm{\Omega }}}^{{\rm{T}}}M{\rm{\Omega }}{{\rm{\Omega }}}^{{\rm{T}}}{\boldsymbol{X}}\\ & = & {{\boldsymbol{\beta }}}^{{\rm{T}}}(-{{\rm{\Omega }}}^{{\rm{T}}}M{\rm{\Omega }}){{\rm{\Omega }}}^{{\rm{T}}}{\boldsymbol{X}}\\ & = & {{\boldsymbol{\beta }}}^{{\rm{T}}}({\rm{\Omega }}M{\rm{\Omega }}){{\rm{\Omega }}}^{{\rm{T}}}{\boldsymbol{X}}\\ & = & {{\boldsymbol{\gamma }}}^{{\rm{T}}}{{\rm{\Omega }}}^{{\rm{T}}}{\boldsymbol{X}},\end{eqnarray*}$$

whence

$$\begin{eqnarray*}&&{\chi }^{\prime }({\boldsymbol{\beta }})={\rm{Tr}}\{\hat{\rho }\mathrm{exp}[{\rm{i}}{\boldsymbol{X}}\cdot {\rm{\Omega }}{\boldsymbol{\gamma }}]\}{\rm{Tr}}\{\hat{\rho }\mathrm{exp}[{\rm{i}}{\boldsymbol{F}}\cdot {\rm{\Omega }}{\boldsymbol{\beta }}]\}.\end{eqnarray*}$$

We now recognise the term on the left as being simply the input characteristic function with its arguments rescaled to ${\boldsymbol{\gamma }}$, viz

$$\begin{eqnarray*}&&\chi ({\boldsymbol{\beta }})\to \chi ({\rm{\Omega }}{M}^{{\rm{T}}}{\rm{\Omega }}{\boldsymbol{\beta }}){\rm{Tr}}\{\hat{\rho }\mathrm{exp}[{\rm{i}}{\boldsymbol{F}}\cdot {\rm{\Omega }}{\boldsymbol{\beta }}]\}.\end{eqnarray*}$$

Since we are already aware of the general form of a Gaussian characteristic function (E3a) we may simply use this formula to find the kernel induced by Gaussian noise, viz

$$\begin{eqnarray*}&&{\chi }_{{\rm{F}}}({\boldsymbol{\beta }})={\rm{Tr}}\{\hat{\rho }\mathrm{exp}[{\rm{i}}{\boldsymbol{F}}\cdot {\rm{\Omega }}{\boldsymbol{\beta }}]\}=\mathrm{exp}\left\{-\displaystyle \frac{1}{2}{{\boldsymbol{\beta }}}^{{\rm{T}}}{{\rm{\Omega }}}^{{\rm{T}}}{V}_{\mathrm{FF}}{\rm{\Omega }}{\boldsymbol{\beta }}\right\},\end{eqnarray*}$$

where ${V}_{\mathrm{FF}}={\mathfrak{R}}\{\langle {\boldsymbol{F}}{{\boldsymbol{F}}}^{{\rm{T}}}\rangle \}$.

Let us now consider a two mode system. If the density matrix is separable (i.e. $\hat{\rho }={\hat{\rho }}_{{\rm{M}}}\otimes {\hat{\rho }}_{{\rm{L}}}$), as we have assumed, the input characteristic function factorises into the form

$$\begin{eqnarray*}&&\chi ({\boldsymbol{\beta }})={\chi }_{{\rm{M}}}({{\boldsymbol{\beta }}}_{{\rm{M}}}){\chi }_{{\rm{L}}}({{\boldsymbol{\beta }}}_{{\rm{L}}}),\end{eqnarray*}$$

where

$$\begin{eqnarray*}&&{\boldsymbol{\beta }}=\left(\begin{array}{c}{{\boldsymbol{\beta }}}_{{\rm{M}}}\\ {{\boldsymbol{\beta }}}_{{\rm{L}}}\end{array}\right).\end{eqnarray*}$$

Let us define the matrices

$$\begin{eqnarray*}{E}_{{\rm{M}}} & = & \left(\begin{array}{cccc}1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\end{array}\right),\\ {E}_{{\rm{L}}} & = & \left(\begin{array}{cccc}0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1\end{array}\right).\end{eqnarray*}$$

With these we may write ${{\boldsymbol{\beta }}}_{{\rm{M}},{\rm{L}}}={E}_{{\rm{M}},{\rm{L}}}{\boldsymbol{\beta }}$. Thus, from the evolution above we can find that

$$\begin{eqnarray*}&&\chi ({\rm{\Omega }}\;{M}^{{\rm{T}}}\;{\rm{\Omega }}{\boldsymbol{\beta }})={\chi }_{{\rm{M}}}({E}_{{\rm{M}}}{\rm{\Omega }}\;{M}^{{\rm{T}}}\;{\rm{\Omega }}{\boldsymbol{\beta }}){\chi }_{{\rm{L}}}({E}_{{\rm{L}}}{\rm{\Omega }}\;{M}^{{\rm{T}}}\;{\rm{\Omega }}{\boldsymbol{\beta }}).\end{eqnarray*}$$

This allows us to write the transformation between input and output characteristic functions as

$$\begin{eqnarray*}&&{\chi }_{{\rm{M}}}({\boldsymbol{\beta }}){\chi }_{{\rm{L}}}({\boldsymbol{\beta }})\to {\chi }_{{\rm{M}}}({E}_{{\rm{M}}}{\rm{\Omega }}\;{M}^{{\rm{T}}}\;{\rm{\Omega }}{\boldsymbol{\beta }}){\chi }_{{\rm{L}}}({E}_{{\rm{L}}}{\rm{\Omega }}\;{M}^{{\rm{T}}}\;{\rm{\Omega }}{\boldsymbol{\beta }}){\chi }_{{\rm{F}}}({\boldsymbol{\beta }}).\end{eqnarray*}$$

The final Wigner function of the mechanical oscillator is obtained by evaluating this at ${{\boldsymbol{\beta }}}^{{\rm{T}}}=({\mathfrak{R}}\{{\beta }_{{\rm{M}}}\},{\mathfrak{I}}\{{\beta }_{{\rm{M}}}\},0,0)$ and applying the appropriate transformation (E1).

Note that this calculation is equivalent to determining the Wigner function by propagating the reduced density matrix of the mechanical mode in the Schrödinger picture.

E.3. Fidelity

The Schumacher fidelity may be recast as

$$\begin{eqnarray*}&&{ \mathcal F }=4\pi \int {{\rm{d}}}^{2}{\boldsymbol{r}}\ {W}^{\prime }({\boldsymbol{r}}){W}_{\mathrm{target}}({\boldsymbol{r}}).\end{eqnarray*}$$