Multimode phononic correlations in a nondegenerate parametric amplifier

The quantum control, detection and manipulation of macroscopic mechanical systems has made enormous strides from its origins in the context of gravitational wave detection [1–3] to current efforts on cavity optomechanics [4–8] and quantum non-demolition measurements [9–11]. With increasing sophistication of experimental techniques and material platforms amenable to such studies, a broader range of questions have come into focus including the use of such mesoscopic mechanical systems for studies of macroscopic entanglement [12–15], out-of-equilibrium thermodynamics [16, 17] and the quantum-to-classical transition [18–20].

A key enabling ingredient for these studies is the realization of mechanical systems with low dissipation and strong, quantum-compatible nonlinear interactions. While mesoscopic mechanical systems exhibit a wide range of mechanical nonlinearities [21], it is typically the case that such nonlinear effects are weak and only arise at large amplitudes of motion or are present in highly dissipative systems. In either scenario, these preclude quantum-limited operation. Alternately, such nonlinear couplings can also be realized through optical mediation [22, 23]. However, the experimental constraints posed by such optically mediated interactions remain challenging to satisfy.

An alternate avenue to combining low intrinsic dissipation and strong nonlinear interactions exploits the notion of reservoir engineering [24–27]—the control of the properties and effective interactions of the system through appropriate design of its environment. Most proposals to date have focused on tuning the properties of an optical reservoir that is coupled to the mechanical system via radiation pressure, a coupling that is typically weak in current optomechanical systems. However, reservoir engineering can also be effected through purely mechanical means, i.e. the interactions between distinct modes of a resonator can be mediated and enhanced by discrete excitations of a massive supporting substrate. This has recently been realized experimentally in an ultrahigh Q membrane resonator in [28], where mechanical parametric nonlinearities, down-conversion and two-mode thermomechanical noise squeezing is demonstrated in a system compatible with established techniques of radiation pressure cooling to the quantum regime [29, 30] and quantum-limited optical detection [31, 32]. Such nonlinear phenomena have also been engineered through a geometric coupling between distinct electromechanical beam resonators [33, 34].

The engineering of strong two-mode parametric nonlinearities in quantum-compatible mechanical resonators [28] has set the stage for the creation of nonclassical mechanical states and the manipulation of phononic fields in a manner akin to quantum optical processes in optical parametric amplifiers and oscillators (OPA/OPOs). In OPA/OPOs, nonclassical correlations between distinct optical fields can be realized at the single photon level as a result of the coherent down-conversion of a high frequency photon into two lower frequency photons. OPA/OPOs have been used to generate squeezed light [35], demonstrate continuous variable EPR entanglement [36, 37] and have several applications in quantum information [38], communication [39, 40] and metrology [41]. Similarly, the realization of quantum-compatible nonlinear phononic processes in mechanical systems offers rich prospects for studies of nonlinear metrology, the robust generation of entangled mechanical states and the quantum dynamics of mesoscopic mechanical systems.

Apart from realizing mechanical analogs of the phenomena in OPA/OPOs, parametric nonlinear process in quantum-compatible mechanical systems also offer unique opportunities in their own right. Unlike in OPA/OPOs where the intracavity fields can only be indirectly detected, mechanical systems allow for the direct nondestructive detection of their displacement—the phononic analog of the intracavity light field. This allows for novel implementations of interferometric schemes in mechanical systems that surpass the standard quantum limit [42]. Furthermore, the extremely low rates of thermalization realized in these mechanical systems [43] open novel avenues for studies of non-equilibrium physics of open quantum systems. For instance, the parametric two-mode nonlinearity described in this work can be harnessed for studies of critical dynamics and entanglement, and the robust generation of non-Gaussian states by appropriate reservoir engineering [44, 45].

Given this range of opportunities, we present a comprehensive description of this two-mode nonlinearity describing various regimes of operation. We emphasize the nature and fidelity of the multimode phononic correlations arising from the nonlinear coupling. While we have focused on nonlinearities that have been engineered through coupling to a reservoir of substrate modes [28], our analysis and the calculated squeezing spectra are also valid for nonlinear interactions engineered through other means [33].

The paper is organized as follows—in section 2 we describe the two-mode nonlinearity arising from the parametric interaction between the resonator and its supporting substrate. This interaction realizes a phononic version of a nondegenerate parametric amplifier involving a substrate excitation (pump) and two resonator modes (signal, idler). Such amplifiers are characterized by a threshold pump amplitude beyond which the system is susceptible to self-oscillation. In section 3 we describe the behavior of this system below threshold. We discuss the onset of two-mode squeezing and calculate the limits of such squeezing in the presence of dissipation, frequency asymmetries and finite pump detuning from parametric resonance. In section 4, we describe the behavior of this system above threshold where the parametric oscillator exhibits amplitude difference squeezing. Finally, in section 5 we describe the crossover regime around the instability threshold and applications of this system to nonlinear metrology and mesoscopic quantum dynamics.

We begin by describing the phenomenology of the two-mode parametric nonlinearity. To motivate the discussion, we consider the physical system described in previous work [28]. The mechanical resonator consists of a silicon nitride (SiN) membrane under high tensile stress that is deposited on a single-crystal silicon substrate. These membranes are a promising optomechanical platform due to their low optical absorption and ultralow dissipation [46, 47]. Their excellent mechanical properties are due to a combination of high intrinsic stress and substrate-induced suppression of anchor loss. This leads to the robust formation of a large number of mechanical modes with low dissipation and high degree of isolation from the environment [43]. Furthermore, the supporting substrate can parametrically mediate multimode interactions within the membrane. As shown in [28], such nonlinear interactions are especially significant when the parametric resonance coincides with a discrete excitation of the substrate, as the coupling strength is now enhanced by the quality factor of the mediating excitation (see figure 1).

Zoom In Zoom Out Reset image size

We describe the parametric two-mode nonlinearity by an interaction of the form $\mathcal{H}=-{{gX}}_{S}{x}_{i}{x}_{j},$ where g parametrizes the strength of the interaction, X_S is the amplitude of the substrate excitation and ${x}_{i,j}$ are the amplitudes of the individual membrane resonator modes.

Within the rotating wave approximation, this results in equations of motion of the form

$$\begin{eqnarray}&&{\ddot{x}}_{i}+{\gamma }_{i}{\stackrel{\dot{}}{x}}_{i}+{\omega }_{i}^{2}{x}_{i}={\displaystyle \frac{1}{{m}_{i}}}({F}_{i}(t)+{\displaystyle \frac{g}{2}}{X}_{S}{x}_{j}^{*})\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\ddot{X}}_{S}+{\gamma }_{S}{\stackrel{\dot{}}{X}}_{S}+{\omega }_{S}^{2}{X}_{S}={\displaystyle \frac{1}{{M}_{S}}}({F}_{S}(t)+{\displaystyle \frac{g}{2}}{x}_{i}{x}_{j})\end{eqnarray} \tag{ 2 }$$

with the corresponding equation for mode j obtained by substituting $i\leftrightarrow j$ in equation (1). Here, we have taken ${x}_{i,j},{X}_{S}$ to denote the (complex) displacements of the individual modes. The external actuating force and thermomechanical noise forces acting on the various modes are together represented by ${F}_{i,j,S}$, and ${\omega }_{i,j,S}$, ${\gamma }_{i,j,S}$ and ${m}_{i,j,S}$ are the eigenfrequencies, mechanical linewidths and the masses of the modes. For now, we assume that ${\omega }_{S}={\omega }_{i}+{\omega }_{j}$, i.e. the pump actuation is at the parametric resonance. Lastly, in keeping with the experimental scenario, we assume that the dissipation rate of the substrate excitation is significantly larger than that of the membrane modes (${Q}_{S}={\omega }_{S}{/\gamma }_{S}\sim {10}^{3}-{10}^{4}$, ${Q}_{i,j}={\omega }_{i,j}{/\gamma }_{i,j}\sim {10}^{7}$).

The coupled equations of motion can be solved using the methods of two timescale perturbation theory [21]. This gives first order coupled equations of the form

$$\begin{eqnarray}&&2{\stackrel{\dot{}}{A}}_{i}={\gamma }_{i}\left[-{A}_{i}+{\rm{i}}{\displaystyle \frac{g}{2}}{\chi }_{i}{A}_{j}^{*}{A}_{S}+{\rm{i}}{\chi }_{i}{\tilde{F}}_{i}(t)\right]\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&2{\stackrel{\dot{}}{A}}_{j}={\gamma }_{j}\left[-{A}_{j}+{\rm{i}}{\displaystyle \frac{g}{2}}{\chi }_{j}{A}_{i}^{*}{A}_{S}+{\rm{i}}{\chi }_{j}{\tilde{F}}_{j}(t)\right]\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&2{\stackrel{\dot{}}{A}}_{S}={\gamma }_{S}\left[-{A}_{S}+{\rm{i}}{\displaystyle \frac{g}{2}}{\chi }_{S}{A}_{i}{A}_{j}+{\rm{i}}{\chi }_{S}\tilde{{F}_{S}}(t)\right],\end{eqnarray} \tag{ 5 }$$

where ${x}_{k}={A}_{k}{\rm{e}}^{-{\rm{i}}{\omega }_{k}t},\;k\in [i,j,S]$. Also, ${\tilde{F}}_{k},\;k\in [i,j,S]$ are the slowly varying (complex) amplitudes of the external forces on the individual modes, and ${\chi }_{k}={({m}_{k}{\omega }_{k}{\gamma }_{k})}^{-1}$ are the magnitudes of the on-resonant susceptibility of the various modes. We have ignored terms such as ${\ddot{A}}_{k},{\gamma }_{i}{\stackrel{\dot{}}{A}}_{k}$ in the slow time approximation.

As the actuating force on the substrate is increased from zero, the substrate (pump) displacement increases in linear proportion until it reaches a threshold amplitude

$$\begin{eqnarray}&&| {A}_{S,\mathrm{cr}}| ={\chi }_{S}{\tilde{F}}_{S,\mathrm{cr}}={\displaystyle \frac{2}{g\sqrt{{\chi }_{i}{\chi }_{j}}}}\propto \sqrt{{\displaystyle \frac{1}{{Q}_{i}{Q}_{j}}}}.\end{eqnarray} \tag{ 6 }$$

Below this critical amplitude, the steady-state amplitudes of the membrane modes (signal, idler), denoted by ${\bar{A}}_{i,j}$, remain at zero and the system realizes a nondegenerate parametric amplifier. Once the parametric drive exceeds the critical value, the membrane modes exhibit an instability to self-oscillation with ${\bar{A}}_{i,j}\ne 0$, and the system realizes a phononic version of an optical parametric oscillator. At the threshold, the system is characterized by a divergent mechanical susceptibility and response time. The behavior of the parametric system in the vicinity of this threshold can be described in terms of a nonequilibrium continuous phase transition.

As can be seen from the above expression, the two-mode coupling is enhanced by the quality factors of the individual resonator modes, resulting in strong multimode interaction strengths even in the presence of low dissipation.

As the pump actuation is further increased, the substrate amplitude remains at the threshold value while the signal and idler amplitudes grow as

$$\begin{eqnarray}&&| {\bar{A}}_{i,j}| ={\displaystyle \frac{2}{g\sqrt{{\chi }_{i,j}{\chi }_{S}}}}\sqrt{\mu -1}.\end{eqnarray} \tag{ 7 }$$

This behavior of the pump, signal and idler mode amplitudes as a function of the normalized drive, $\mu =\frac{| {F}_{S}| }{| {F}_{S,\mathrm{cr}}| }$ is shown in figure 2.

Zoom In Zoom Out Reset image size

This parametric process can be viewed as the down-conversion of one phonon from the pump mode to a pair of phonons, one in each of the signal and idler modes. At the quantum limit, this results in the entanglement of motion of these modes. In the classical regime, this down-conversion is manifest as correlations in the mechanical motion of the two modes. Below threshold, the nonlinear interaction realizes a parametric amplifier with a gain that is dependent on the relation between the phases of the resonator modes and that of the pump. This phase-dependent gain results in thermomechanical squeezing of composite quadratures formed from linear combinations of the quadratures of the individual mechanical modes [28, 33].

Above threshold, the rate of phonon down-conversion exceeds the intrinsic loss from either resonator mode, leading to self-oscillation. In this regime, the correlated production of down-converted phonons manifests as a reduction in the variance (squeezing) of the difference in the amplitude fluctuations of the signal and idler modes. This is the thermomechanical analog of intensity difference squeezing in optical parametric oscillators.

In this section, we evaluate the dynamics of the nondegenerate parametric amplifier that arise for pump actuation below the instability threshold. In this regime, the system exhibits correlations between the displacements of the signal and idler modes. These correlations are manifest as two-mode squeezing of a combined quadrature composed from the individual modes.

In general, the correlations can be obtained using the coupled equations for the membrane modes in the simultaneous presence of a classical pump actuation and thermomechanical Langevin forces. The detailed derivation of the noise spectra is given in appendix A and we only briefly outline the procedure here. The analysis of the thermal fluctuation spectra is analogous to that for quantum fluctuations of optical parametric amplifiers [48–50]. This similarity is due to both quantum and thermal fluctuations being approximated as Gaussian. In addition, both the thermal Langevin forces in mechanical systems and the quantum noise in OPA/OPOs are typically assumed to be Markovian. This approximation has been shown to breakdown in the thermal baths of mechanical resonators [44]. The consequences of this breakdown and the nature of the squeezing spectra in the presence of non-Markovian corrections to thermal fluctuations will be described elsewhere.

We separate the mean displacement and fluctuations about the mean, by writing ${x}_{i,j}=({\bar{A}}_{i,j}+{\scr{\Delta}} {A}_{i,j}){\rm{e}}^{-{\rm{i}}{\omega }_{i,j}t}$ with $\langle {\scr{\Delta}} {A}_{i,j}\rangle =0$. Further, the fluctuations are decomposed into their quadrature components as ${\scr{\Delta}} {A}_{i,j}={\scr{\Delta}} {\alpha }_{i,j}+{\rm{i}}{\scr{\Delta}} {\beta }_{i,j}$. Below threshold, the fluctuations of the individual quadratures are given by equations (A.4), (A.5) of appendix A with ${\bar{A}}_{i,j}=0$.

We then define cross-quadratures constructed from $\{{\alpha }_{i,j},{\beta }_{i,j}\}$, here normalized to their respective thermomechanical amplitudes, according to the relations

$$\begin{eqnarray}&&{x}_{\pm }=({\alpha }_{i}\pm {\alpha }_{j})/\sqrt{2}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{y}_{\pm }=({\beta }_{i}\pm {\beta }_{j})/\sqrt{2}.\end{eqnarray} \tag{ 9 }$$

The two-mode correlations are manifest as amplification and squeezing of the above quadratures. We represent the fluctuations in these cross-quadratures, along with the fluctuation of the substrate mode in the form of the column vectors ${\bf{X}}={\left({\scr{\Delta}} {x}_{+},{\scr{\Delta}} {x}_{-},{\scr{\Delta}} {x}_{S}\right)}^{T}$, ${\bf{Y}}={\left({\scr{\Delta}} {y}_{+},{\scr{\Delta}} {y}_{-},{\scr{\Delta}} {y}_{S}\right)}^{T}$. These are related to the original quadrature fluctuations ${\scr{\Delta}} \mathop{\alpha }\limits^{\rightarrow }$, ${\scr{\Delta}} \mathop{\beta }\limits^{\rightarrow }$ via

$$\begin{eqnarray}{\bf{X}}={\bf{R}}{\scr{\Delta}} \mathop{\alpha }\limits^{\rightarrow };\quad {\bf{Y}}={\bf{R}}{\scr{\Delta}} \mathop{\beta }\limits^{\rightarrow };\quad {\bf{R}}={\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{ccc}1 & 1 & 0\\ 1 & -1 & 0\\ 0 & 0 & 1\end{array}\right).\end{eqnarray} \tag{ 10 }$$

Finally, the correlations of the cross quadratures are obtained through a corresponding transformation of the spectral densities ${{\bf{S}}}_{\alpha /\beta }(\omega )$, for example

$$\begin{eqnarray}&&{{\bf{S}}}_{X}(\omega )=\left\langle{\bf{X}}(\omega ){\bf{X}}{(\omega )}^{\dagger }\right\rangle={\bf{R}}{{\bf{S}}}_{\alpha }{{\bf{R}}}^{T}\end{eqnarray} \tag{ 11 }$$

along with the analogous equation for ${\bf{Y}}$ ($X\to Y\quad \&\quad \alpha \to \beta$ in the above). As shown in appendix A, the degree of squeezing is obtained through the variances of these cross quadratures by integrating the spectra in equation (11) over the measurement bandwidth.

While this outlines the basis of the calculation, the appearance of such two-mode correlations can be intuitively seen as the result of a coherent interference between the response of the individual resonator modes due to thermomechanical noise and the down-converted field arising from the two-mode nonlinearity. This interference results in a reduction of thermomechanical motion along one quadrature at the expense of amplified motion in the orthogonal quadrature.

In general, the degree of two-mode noise squeezing is sensitive to various experimental considerations such as the 'loss asymmetry' arising from mismatched dissipation rates of the individual membrane modes, their frequency difference and the detuning of the parametric drive from the two-mode resonance. Below, we evaluate the effect of these considerations on the degree of squeezing and find the robust formation of two-mode squeezed states for a wide range of experimental parameters.

3.1. Matched losses and frequencies

For notational convenience, we introduce the loss asymmetry parameter ${{\scr{\Delta}} }_{\gamma }=({\gamma }_{i}-{\gamma }_{j})/({\gamma }_{i}+{\gamma }_{j})$ and frequency mismatch parameter ${{\scr{\Delta}} }_{\omega }=({\omega }_{i}-{\omega }_{j})/({\omega }_{i}+{\omega }_{j})$. To build an intuition for quadrature squeezing below threshold, we first consider the simplest case of distinct resonator modes with identical frequencies (${\omega }_{i}={\omega }_{j}=\omega$ or ${{\scr{\Delta}} }_{\omega }=0$) and identical dissipation rates (${\gamma }_{i}={\gamma }_{j}=\gamma$ or ${{\scr{\Delta}} }_{\gamma }=0$).

For this case, the evolution matrices in equation (A.6) of appendix A reduce to

$$\begin{eqnarray}{{\bf{M}}}_{\alpha /\beta }={\displaystyle \frac{1}{2}}\left(\begin{array}{ccc}-\gamma & \mp \gamma \mu & 0\\ \mp \gamma \mu & -\gamma & 0\\ 0 & 0 & -{\gamma }_{S}\end{array}\right).\end{eqnarray} \tag{ 12 }$$

The spectral density of fluctuations of the collective quadratures is evaluated from equations (A.7) and (11) to yield

$$\begin{eqnarray}{{\bf{S}}}_{{\bf{X}}/Y}={\displaystyle \frac{2}{\pi }}\left(\begin{array}{ccc}\frac{\gamma }{\left({\gamma }^{2}{(1\pm \mu )}^{2}+4{\omega }^{2}\right)} & 0 & 0\\ 0 & \frac{\gamma }{\left({\gamma }^{2}{(1\mp \mu )}^{2}+4{\omega }^{2}\right)} & 0\\ 0 & 0 & \frac{{\gamma }_{S}}{\left({\gamma }_{S}^{2}+4{\omega }^{2}\right)}\end{array}\right).\end{eqnarray} \tag{ 13 }$$

The variances of the normalized collective quadratures of the mechanical modes are then given by

$$\begin{eqnarray}&&{\sigma }_{{x}_{\pm },{x}_{\pm }}={\displaystyle \frac{1}{1\pm \mu }}={\sigma }_{{y}_{\mp },{y}_{\mp }}.\end{eqnarray} \tag{ 14 }$$

We see that ${x}_{-},{y}_{+}$ are amplified quadratures with variances that grow as $\mu \to 1$, while ${x}_{+},{y}_{-}$ are squeezed quadratures showing reduction in the variance below the thermomechanical limit. For pump actuation close to parametric threshold ($\mu \to 1$), we obtain a peak noise squeezing of $\frac{1}{2}$, similar to the bound observed in optical parametric amplifiers [51]. The degree of two-mode squeezing versus parametric drive is shown in figure 3.

Zoom In Zoom Out Reset image size

The existence of a bound for the peak noise squeezing can be intuited by looking at the equations of motion for the collective quadratures, through a rotation of the original equations of motion, i.e.

$$\begin{eqnarray}&&\stackrel{\dot{}}{{\bf{X}}}={{\bf{M}}}_{X}{\bf{X}}+{\bf{v}}{\prime }_{X}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\stackrel{\dot{}}{{\bf{Y}}}={{\bf{M}}}_{Y}{\bf{Y}}+{\bf{v}}{\prime }_{Y},\end{eqnarray} \tag{ 16 }$$

where

$$\begin{eqnarray*}{{\bf{M}}}_{X/Y}={\bf{R}}{{\bf{M}}}_{\alpha /\beta }{{\bf{R}}}^{T}=-{\displaystyle \frac{1}{2}}\left(\begin{array}{ccc}\gamma (1\pm \mu ) & 0 & 0\\ 0 & \gamma (1\mp \mu ) & 0\\ 0 & 0 & {\gamma }_{s}\end{array}\right)\end{eqnarray*}$$

and ${\bf{v}}\prime ={\bf{Rv}}$ represents thermomechanical noise forces.

Two-mode squeezing arises from the fact that while the thermomechanical noise forces remain the same in the presence or absence of the parametric drive, the dissipation of the squeezed quadrature goes to twice its bare value near parametric threshold. This results in a reduction in the variance of the squeezed quadrature by a factor of 2. Simultaneously, the decay rate of the amplified quadrature goes to zero, signaling the onset of the parametric instability. The onset of the instability thus explains the 3 dB bound for squeezing through this parametric process.

The variances of the cross-quadratures are calculated by integrating the noise over all frequencies, i.e. they are measured with infinite bandwidth. As the measurement bandwidth decreases, the peak noise squeezing approaches 6 dB, i.e.

$$\begin{eqnarray}&&{\displaystyle \frac{{S}_{{x}_{+},{x}_{+}}(\omega =0,\mu )}{{S}_{{x}_{+},{x}_{+}}(\omega =0,\mu =0)}}={\displaystyle \frac{1}{{(1+\mu )}^{2}}}\stackrel{\mu \to 1}{\;\rightarrow \;}{\displaystyle \frac{1}{4}}.\end{eqnarray} \tag{ 17 }$$

3.2. Effect of mismatched frequencies and loss rates

We now consider the case where the frequencies and the damping rate of the two resonators are not the same, i.e. ${{\scr{\Delta}} }_{\gamma },{{\scr{\Delta}} }_{\omega }\ne 0$. In this situation, the collective quadratures ${x}_{\pm }$ and ${y}_{\pm }$ are no longer decoupled from each other. This coupling between the collective quadratures has been noted to be detrimental to entanglement [52], and backaction evasion [23] protocols. It also results in a degradation of peak two-mode thermomechanical squeezing. This is a result of the coupled quadratures ${x}_{+},{x}_{-}$ and ${y}_{-},{y}_{+}$ being respectively amplified and squeezed in the presence of the parametric drive.

As before, the variances of the collective quadratures are obtained using equation (A.7), subsequent rotation using equation (11) and integration over frequency. These are now given by

$$\begin{eqnarray}\begin{array}{l}{\sigma }_{{y}_{\pm },{y}_{\pm }} & = & {\sigma }_{{x}_{\mp },{x}_{\mp }}\\ & = & {\displaystyle \frac{1}{1-{\mu }^{2}}}\left\{1+{\mu }^{2}{\displaystyle \frac{{{\scr{\Delta}} }_{\omega }({{\scr{\Delta}} }_{\omega }-{{\scr{\Delta}} }_{\gamma })}{1-{{\scr{\Delta}} }_{\omega }^{2}}}\pm \mu \sqrt{{\displaystyle \frac{1-{{\scr{\Delta}} }_{\gamma }^{2}}{1-{{\scr{\Delta}} }_{\omega }^{2}}}}\right\}\end{array}\end{eqnarray} \tag{ 18 }$$

with the cross correlation between (${x}_{+},{x}_{-}$), (${y}_{+},{y}_{-}$) given by,

$$\begin{eqnarray}&&{\sigma }_{{y}_{+},{y}_{-}}={\sigma }_{{x}_{+},{x}_{-}}={\displaystyle \frac{{\mu }^{2}\left({{\scr{\Delta}} }_{\omega }-{{\scr{\Delta}} }_{\gamma }\right)}{2\left(1-{\mu }^{2}\right)\left(1-{{\scr{\Delta}} }_{\omega }^{2}\right)}}.\end{eqnarray} \tag{ 19 }$$

The variances of the amplified and squeezed collective quadratures for the case where (${{\scr{\Delta}} }_{\gamma }\ne 0,{{\scr{\Delta}} }_{\omega }\ne 0$, ${{\scr{\Delta}} }_{\gamma }\ne {{\scr{\Delta}} }_{\omega }$), are shown in figure 3. As can be seen, the presence of loss asymmetry or a frequency mismatch results in a degradation of noise squeezing. In this case, optimal squeezing obtained for a parametric drive that is significantly below the instability threshold. We also find that the coupling between the amplified and squeezed quadratures leads to a divergence of the squeezed quadrature at the instability threshold.

The dependence of the peak squeezing on the loss asymmetry and frequency mismatch parameters are summarized in figure 4. Figure 4(a) shows a plot of the peak squeezing as a function of the loss asymmetry for the case of distinct mechanical modes with the same frequency (${{\scr{\Delta}} }_{\omega }=0$), showing a linear degradation of peak squeezing with loss asymmetry. Correspondingly, figure 4(b) shows the peak squeezing as a function of the frequency mismatch parameter (${{\scr{\Delta}} }_{\omega }$) for the case of no loss asymmetry (${{\scr{\Delta}} }_{\gamma }=0$).

Zoom In Zoom Out Reset image size

Importantly, as can be seen in figure 4(c), we find that the 3 dB squeezing bound can be regained even in the presence of loss asymmetry as long as ${{\scr{\Delta}} }_{\gamma }={{\scr{\Delta}} }_{\omega }$. In this case, the normalized cross correlations between amplified and squeezed collective quadratures vanish, and equation (18) reduces to equation (14). This results in a noise squeezing that is identical to that for the case of symmetric losses and frequencies.

3.3. Effect of pump detuning

Finally, we consider the thermomechanical squeezing bound when the pump drive is detuned from parametric resonance. In this case, the drive frequency is given by ${\omega }_{d}={\omega }_{S}+{\scr{\Delta}}$, where Δ is the drive detuning. This case is of interest since a detuned parametric drive introduces a dynamic coupling between the ${\scr{\Delta}} \mathop{\alpha }\limits^{\rightarrow }$ and ${\scr{\Delta}} \mathop{\beta }\limits^{\rightarrow }$ quadratures. In turn, this leads to correlations between the amplified and squeezed quadratures. Due to these correlations, the amplified quadrature contains information about the squeezed quadrature that can be used for enhanced localization through weak measurements and optimal estimation [53]. Thus, a detuned parametric drive can allow for enhanced noise squeezing in the presence of feedback. An additional point of interest in this case is that special choices of the drive detuning lead to some of the collective quadratures becoming quantum non-demolition observables [52].

The equations satisfied by the slowly varying complex amplitudes (${A}_{i,j,S}$) are again given by equations (3)–(5), with the only difference now being that the pump actuation ${\tilde{F}}_{S}(t)$ is a slowly varying function of time, ${\tilde{F}}_{S}(t)=| {\tilde{F}}_{S}| {\rm{e}}^{-{\rm{i}}{\scr{\Delta}} t}$. The pump amplitude resulting from this drive force, ${\bar{A}}_{S}(t)$ is given by

$$\begin{eqnarray}&&{\bar{A}}_{S}={\rm{i}}{\chi }_{S}{\tilde{F}}_{S}(t)={\rm{i}}{\chi }_{S}| {\tilde{F}}_{S}| {\rm{e}}^{-{\rm{i}}{\scr{\Delta}} t}={\rm{i}}| {\bar{A}}_{S}| {\rm{e}}^{-{\rm{i}}{\scr{\Delta}} t}.\end{eqnarray} \tag{ 20 }$$

Here, as we are most interested in pump detunings that are comparable to the resonator linewidths, we have made the assumption that ${\scr{\Delta}} \ll {\gamma }_{S}$ and that the pump amplitude is hence related to the instantaneous parametric actuation through the on-resonant susceptibility.

The two-mode correlations are computed in appendix B, and are given for ${{\scr{\Delta}} }_{\gamma }={{\scr{\Delta}} }_{\omega }=0$ by

$$\begin{eqnarray}\begin{array}{l}{\sigma }_{{x}_{\pm },{x}_{\pm }} & = & \left({\displaystyle \frac{{k}_{B}T\gamma }{m{\omega }^{2}}}\right)\left[{\displaystyle \frac{\left({{\scr{\Delta}} }^{2}+{\gamma }^{2}{(1\mp \mu )}^{2}-{\lambda }_{-}^{2}\right)}{{\lambda }_{+}{\lambda }_{-}\left({\lambda }_{+}+{\lambda }_{-}\right)}}+{\displaystyle \frac{1}{{\lambda }_{+}}}\right]\\ & = & {\sigma }_{{y}_{\mp },{y}_{\mp }}\end{array}\end{eqnarray} \tag{ 21 }$$

and

$$\begin{eqnarray}&&{\sigma }_{{x}_{\pm },{y}_{\pm }}=\left({\displaystyle \frac{{k}_{B}T\gamma }{m{\omega }^{2}}}\right){\displaystyle \frac{2{\scr{\Delta}} \gamma \mu }{{\lambda }_{+}{\lambda }_{-}\left({\lambda }_{+}+{\lambda }_{-}\right)}}={\sigma }_{{y}_{\pm },{x}_{\pm }},\end{eqnarray} \tag{ 22 }$$

where ${\lambda }_{\pm }^{2}={\gamma }^{2}\left(1+{\mu }^{2}\right)-{{\scr{\Delta}} }^{2}\pm 2\gamma \sqrt{{\gamma }^{2}{\mu }^{2}-{{\scr{\Delta}} }^{2}}$.

Note that ${x}_{+}$ and ${y}_{-}$ are amplified quadratures while ${x}_{-}$ and ${y}_{+}$ are squeezed. As mentioned earlier, a nonzero detuning introduces correlations between (${x}_{+},{y}_{+}$) and (${x}_{-},{y}_{-}$), i.e. between amplified and squeezed quadratures. This is distinct from the case of loss asymmetry or frequency mismatch considered previously, where correlations were introduced between (${x}_{+},{x}_{-}$) and (${y}_{+},{y}_{-}$). As can be seen from the above expressions, these correlations between the ${x}_{\pm }$ and ${y}_{\pm }$ quadratures are proportional to the drive detuning.

This coupling between amplified and squeezed quadratures also results in a decrease in the peak squeezing at non-zero detunings, as can be seen in figure 5(a). As expected, for the detuned case, the amplified quadrature diverges at $\mu =\sqrt{1+{({\scr{\Delta}} /\gamma )}^{2}}$, the instability threshold for ${\scr{\Delta}} \ne 0$.

Zoom In Zoom Out Reset image size

The peak squeezing as a function of the detuning, normalized with respect to the decay rate (for no loss asymmetry or frequency mismatch) is shown in figure 5(b). We see that squeezing is almost completely lost when the detuning becomes comparable to the linewidth of the signal and idler modes.

For pump actuation above the parametric instability threshold, the two-mode nonlinearity results in parametric self-oscillation of the individual membrane modes. In this regime, the correlated production of down-converted phonons in the signal and idler modes results in a reduction of fluctuations in the amplitude difference between these two modes. This is the phononic version of intensity difference squeezing observed in optical parametric oscillators [49]. As in the optical case, this form of difference squeezing is of interest in nonlinear interferometric schemes [42] capable of measurement sensitivities surpassing the standard quantum limit.

In this section, we discuss the dynamics of the two-mode nonlinearity above threshold and compute the amplitude fluctuations of the mechanical modes. Additionally, we discuss the fluctuations of their phases. While energy conservation dictates that the sum of the frequencies of the signal and idler modes equal that of the pump drive, it is also the case that the sum of the phases of the signal and idler modes above threshold remain locked to the phase of the pump drive. This follows from equation (5), which in steady state can be rewritten in terms of the normalized drive μ and the drive force phase ${\phi }_{S}$ (${F}_{S}=| {F}_{S}| {\rm{e}}^{{\rm{i}}{\phi }_{S}}$) as

$$\begin{eqnarray}&&-(\mu -1){\rm{e}}^{{\rm{i}}{\phi }_{S}}={\displaystyle \frac{g}{2}}{\displaystyle \frac{{\chi }_{S}}{| {A}_{S,\mathrm{cr}}| }}{A}_{i}{A}_{j}.\end{eqnarray} \tag{ 23 }$$

Defining the phases of the resonator modes as ${A}_{i,j}={\rm{i}}| {A}_{i,j}| {\rm{e}}^{{\rm{i}}{\phi }_{i,j}}$, we see that ${\phi }_{i}+{\phi }_{j}={\phi }_{S}$. However, the difference in the phases is not constrained and is free to fluctuate.

As before, we quantify these fluctuations in the amplitudes and phases about their steady state values by decomposing the complex amplitude fluctuations ${\scr{\Delta}} \mathop{A }\limits^{\rightarrow }$ into ${\scr{\Delta}} {\alpha }_{i}$ and ${\scr{\Delta}} {\beta }_{i}$ quadratures. The equations of motion of ${\scr{\Delta}} \mathop{\alpha }\limits^{\rightarrow }$ and ${\scr{\Delta}} \mathop{\beta }\limits^{\rightarrow }$ are given by equations (A.4), (A.5) in appendix A after substituting $| {\bar{A}}_{S}| =\sqrt{\frac{{\gamma }_{i}{\gamma }_{j}}{{\kappa }_{i}{\kappa }_{j}}}$ and ${\bar{A}}_{i,j}=\sqrt{\frac{{\gamma }_{i,j}{\gamma }_{S}}{{\kappa }_{i,j}{\kappa }_{S}}}\sqrt{\mu -1}$. Here, we have defined the coupling parameter ${\kappa }_{k}=\frac{g}{2{m}_{k}{\omega }_{k}}$.

We choose the drive and resonator mode phases such that ${\bar{\phi }}_{i,j}={\phi }_{S}=0$. Given this choice of phases, the complex mean amplitudes ${\bar{A}}_{i,j}$ and the fluctuations are as shown in the inset of figure 6. With this convention, the amplitude fluctuations are given by ${\scr{\Delta}} {\beta }_{i,j}$. The fluctuations in the phase are obtained from ${\scr{\Delta}} {\alpha }_{i,j}$ through ${\scr{\Delta}} {\phi }_{i,j}=\frac{{\scr{\Delta}} {\alpha }_{i,j}}{{\bar{A}}_{i,j}}$.

Zoom In Zoom Out Reset image size

Similar to the below threshold case, the correlations of the fluctuations in the signal and idler modes resulting from the pump drive are manifest in combined quadratures

$$\begin{eqnarray}&&{\scr{\Delta}} {x}_{\pm }={\displaystyle \frac{1}{\sqrt{2}}}({\scr{\Delta}} {\alpha }_{i}\pm {\scr{\Delta}} {\alpha }_{j})\propto {\scr{\Delta}} {\phi }_{\pm }\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{\scr{\Delta}} {y}_{\pm }={\displaystyle \frac{1}{\sqrt{2}}}({\scr{\Delta}} {\beta }_{i}\pm {\scr{\Delta}} {\beta }_{j})={\scr{\Delta}} {R}_{\pm },\end{eqnarray} \tag{ 25 }$$

where ${\scr{\Delta}} {R}_{\pm }$ are the amplitude sum and difference quadratures and ${\scr{\Delta}} {\phi }_{\pm }$ are the phase sum and difference quadratures. The spectrum of fluctuations of ${\scr{\Delta}} {\alpha }_{\pm }$ and ${\scr{\Delta}} {\beta }_{\pm }$ are given by equation (A.7).

In contrast to the dynamics below threshold, the fluctuations of the pump mode above threshold have an influence on the correlations generated between the membrane modes. On the one hand, the substrate fluctuations are smaller than those of the membrane modes by a term proportional to the respective mass ratios. Thus, one might assume that these tiny fluctuations should have a negligible influence on the membrane modes. However, the substrate mode fluctuations affect the membrane modes through terms that are proportional to the steady state amplitudes of the latter. These amplitudes are larger than the pump mode amplitude by the mass ratio of the substrate and membrane modes. These ratios cancel each other and lead to pump mode thermal fluctuations influencing the signal and idler modes through terms that are of the same order as the coupling between the signal and idler modes. Thus, the degree of amplitude difference squeezing is independent of the resonator to substrate mass ratio.

Consistent with the experimental system under consideration [28], we assume that the damping rate of the substrate excitation is 3–4 orders of magnitude larger than those of the membrane. In this regime, the pump fluctuations respond instantaneously to those of the membrane modes and can thus be adiabatically eliminated. We use this to simplify the analysis and ignore the time derivative of the pump fluctuations (${\scr{\Delta}} {\stackrel{\dot{}}{A}}_{S}$) in equation (A.1). Aside from this modification, we extract the fluctuations of the signal and idler modes as before.

4.1. Matched losses and frequencies

We first consider the case where the damping rates and the frequencies are matched. In this limit, we obtain the following spectral densities for the collective quadratures normalized with respect to the thermal motion amplitude

$$\begin{eqnarray}{{\bf{S}}}_{Y}(\omega )={\displaystyle \frac{1}{2\pi }}\left(\begin{array}{cc}\frac{\gamma \mu }{\left({\gamma }^{2}{(-1+\mu )}^{2}+{\omega }^{2}\right)} & 0\\ 0 & \frac{\gamma }{\left({\gamma }^{2}+{\omega }^{2}\right)}\end{array}\right)\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}{{\bf{S}}}_{X}(\omega )={\displaystyle \frac{1}{2\pi }}\left(\begin{array}{cc}\frac{\gamma \mu }{\left({\gamma }^{2}{\mu }^{2}+{\omega }^{2}\right)} & 0\\ 0 & \frac{\gamma }{{\omega }^{2}}\end{array}\right).\end{eqnarray} \tag{ 27 }$$

We obtain the variances of the fluctuations by integrating the spectra. For the Y quadrature, which relates to amplitude fluctuations, these evaluate to

$$\begin{eqnarray}&&{\sigma }_{{y}_{+},{y}_{+}}={\displaystyle \frac{\mu }{2(\mu -1)}}\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{\sigma }_{{y}_{-},{y}_{-}}={\displaystyle \frac{1}{2}}.\end{eqnarray} \tag{ 29 }$$

These variances are plotted as a function of the parametric drive in figure 6.

Above threshold, the amplitude difference between the signal and idler modes is always half the thermal variance. This is the mechanical analogue of intensity difference squeezing seen in optical parametric oscillators. We find that while the individual amplitudes are sensitive to fluctuations of the pump mode, the amplitude difference is insensitive to fluctuations of the pump mode and the degree of squeezing is independent of the pump drive. On the other hand, the variance of the amplitude sum diverges as $\mu \to {1}^{+}$ and decreases with increasing drive, approaching half the thermal variance as $\mu \to \infty$.

The other notable feature of above threshold dynamics is that the phase difference between the signal and idler modes is unspecified and hence free to fluctuate. The fluctuation in the phase difference is given by

$$\begin{eqnarray}&&{S}_{{\phi }_{-},{\phi }_{-}}(\omega )=\left({\displaystyle \frac{{x}_{\mathrm{th}}^{2}}{A{(\mu )}^{2}}}\right){S}_{{x}_{-},{x}_{-}}(\omega )\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&=\;\left({\displaystyle \frac{{x}_{\mathrm{th}}^{2}}{A{(\mu )}^{2}}}\right)\left({\displaystyle \frac{\gamma }{2\pi {\omega }^{2}}}\right),\end{eqnarray} \tag{ 31 }$$

where $A(\mu )$ is the amplitude of the membrane modes (identical for the case of matched loss rates), ${x}_{\mathrm{th}}^{2}=\frac{{k}_{B}T}{m{\omega }^{2}}$ is the thermal variance. The integral of the fluctuation spectrum diverges as ${\omega }^{-2}$, indicating that the difference phase undergoes diffusion. The time scale τ for this diffusion can be estimated by calculating the variance while imposing a low frequency cut off ($\frac{2\pi }{\tau }$) to the integral of the spectral density, i.e.

$$\begin{eqnarray}&&\langle {\scr{\Delta}} {\phi }_{-}^{2}\rangle =2{{\displaystyle \int }}_{{\displaystyle \frac{2\pi }{\tau }}}^{\infty }{S}_{{\phi }_{-},{\phi }_{-}}(\omega ){\rm{d}}\omega \end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&=\;\left({\displaystyle \frac{{x}_{\mathrm{th}}^{2}}{A{(\mu )}^{2}}}\right){\displaystyle \frac{\gamma \tau }{2{\pi }^{2}}}.\end{eqnarray} \tag{ 33 }$$

If the parametric nonlinearity were used to actuate the membrane modes to an amplitude of ${10}^{3}\times {({x}_{\mathrm{th}}^{2})}^{1/2}$, a phase fluctuation of 1 mrad would require a measurement duration of 10 ringdown periods or around 100 s for the modes considered in [28]. Thus, while ${S}_{{\phi }_{-},{\phi }_{-}}(\omega )$ diverges, it does not necessarily lead to large fluctuations of the difference phase over experimental time scales.

4.2. Effect of mismatched frequencies and loss rates

Finally, we discuss the dynamics of the two-mode nonlinearity in the vicinity of the parametric instability. The crossover regime is most conveniently portrayed by evaluating the variance of ${y}_{\pm }$ quadratures in the regimes below and above threshold. In the former regime, these quadratures represent the two-mode correlations arising from phase-sensitive parametric amplification. In the latter regime, these correspond to the sum and difference amplitude fluctuations of the two membrane modes. These are shown in figure 7. In the general case of mismatched dissipation rates (${{\scr{\Delta}} }_{\gamma }\ne 0$), both these quadratures exhibit diverging variances at the instability.

Zoom In Zoom Out Reset image size

As for the divergent phase diffusion discussed earlier, finite time effects need to be considered to interpret the divergent steady-state variances depicted in figure 7. For the low frequency, ultrahigh quality factor resonators considered in this work, the divergent response time in the vicinity of the instability threshold can result in inordinately long thermalization times ($\sim {10}^{4}-{10}^{5}\;{\rm{s}}$). For typical measurement durations (∼100 s), the measured squeezing spectra can deviate appreciably from the spectra computed in steady state. Expectedly, these deviations are most significant for parametric actuation around the instability threshold ($\mu \sim 1.0$).

The variances measured over a finite time ${\tau }_{m}$ are extracted by truncating the integral of the relevant spectral densities by the time of measurement, i.e.

$$\begin{eqnarray}&&{\boldsymbol{\sigma }}_{\alpha ,\beta }=2{{\displaystyle \int }}_{{\displaystyle \frac{2\pi }{{\tau }_{m}}}}^{\infty }{{\bf{S}}}_{\alpha ,\beta }(\omega ){\rm{d}}\omega .\end{eqnarray} \tag{ 34 }$$

These variances, computed for the parameters in [28], result in the solid black curves shown in figure 8. We see that the singularities in the amplified and squeezed quadratures seen in the steady state variances (solid blue and red lines) are washed out at finite measurement durations.

Zoom In Zoom Out Reset image size

The regime around $\mu =1$ is of interest since the system exhibits mechanical bistability and a hysteretic response due to the divergent mechanical susceptibility and diverging response time. This regime offers a clean, mesoscopic and mechanical realization of a second-order phase transition [54–57] and out-of-equilibrium quantum dynamics on experimentally accessible timescales. Furthermore, quantum tunneling between bistable mechanical states has been discussed in the context of a Duffing nonlinearity as a means of accessing a quantum-to-classical transition [20]. In the system considered here and in [28], the presence of a strong two-mode nonlinearity and compatibility with optomechanical cooling together imply that similar effects can be accessed even in the regime of low phonon number and motion on the scale of mechanical zero-point fluctuations. The nature of this nonequilibrium phase transition and the response to critical fluctuations as the phonon occupancy is gradually reduced by optomechanical cooling will be described elsewhere.

We also note the close correspondence between the squeezing spectra arising from the two-mode nonlinearity and the properties of the reservoir to which the resonator is coupled. In particular, while we have considered the squeezing spectra in the presence of a Markovian (thermal) reservoir in this work, it is also known that intrinsic defects or two-level systems (TLS) in amorphous resonators such as SiN can couple to mechanical motion [58]. Furthermore, a reservoir of such TLS can acquire non-Markovian properties in certain regimes. In addition to studying the effect of such non-Markovian fluctuations on the second order phase transition and the two-mode correlations near the instability, it is an intriguing prospect to use the exquisite sensitivity of this two-mode nonlinearity as an amplifier of such reservoir interactions to shed light on the intrinsic material properties of the membrane resonator.

In summary, we describe a phononic nondegenerate parametric amplifier that is realized in a membrane resonator through a substrate-mediated nonlinearity. Motivated by recent work [28, 33], we discuss the creation of multimode phononic correlations arising from this parametric interaction and compute the squeezing spectra in the presence of thermomechanical noise. We address various points of experimental relevance including the presence of frequency mismatches and asymmetric dissipation rates between the resonator modes. We find the robust presence of two-mode phononic correlations for a wide range of experimental parameters.

Below the threshold for parametric instability, this system exhibits two-mode noise squeezing of collective quadratures composed from the individual resonator modes. This regime is conducive to back-action evading schemes for quantum-enhanced metrology. Above threshold, the system exhibits amplitude difference squeezing due to the correlated production of down-converted phonons in both resonator modes. This regime is promising for nonlinear measurement schemes capable of measurement sensitivities surpassing the standard quantum limit.

In the crossover regime between these two limits, the response of the system near the parametric instability threshold is characterized by a divergent mechanical susceptibility and a diverging response time. In this regime, the system exhibits mechanical bistability, a hysteretic response and critical slowing down. This regime is of interest as it offers a mechanical realization of a second order phase transition for the investigation of nonequilibrium critical dynamics, the quantum-to-classical transition and the influence of a non-Markovian reservoir on such critical phenomena. These results will be discussed elsewhere. Due to the combined presence of low intrinsic dissipation and optomechanical compatibility, such nonlinear mechanical systems are promising for extending these studies from the classical realm deep into the quantum regime.

This work was supported by the DARPA QuASAR program through a grant from the ARO, the ARO MURI on non-equilibrium Many-body Dynamics (63834-PH-MUR), an NSF INSPIRE award and the Cornell Center for Materials Research with funding from the NSF MRSEC program (DMR-1120296). We acknowledge valuable discussions with H F H Cheung. M V acknowledges support from the Alfred P Sloan Foundation.

For both the below and the above threshold regimes, we obtain the correlations that develop between the resonator modes in the presence of the pump drive, by analyzing the coupled equations for the resonator modes under the influence of a classical actuation of the pump/substrate mode along with thermomechanical Langevin noise forces acting on the membrane and substrate modes.

We distinguish between the mean displacement and the fluctuations by writing ${x}_{i,j}=({\bar{A}}_{i,j}+{\scr{\Delta}} {A}_{i,j}){\rm{e}}^{-{\rm{i}}{\omega }_{i,j}t}$ where $\langle {\scr{\Delta}} {A}_{i,j}\rangle =0$. The coupled equations for the fluctuations can be written as

$$\begin{eqnarray*}&&\begin{array}{l}2\left(\begin{array}{c}{\scr{\Delta}} {\stackrel{\dot{}}{A}}_{i}\\ {\scr{\Delta}} {\stackrel{\dot{}}{A}}_{j}\\ {\scr{\Delta}} {\stackrel{\dot{}}{A}}_{S}\end{array}\right) & = & \left(\begin{array}{ccc}-{\gamma }_{i} & 0 & {\rm{i}}{\kappa }_{i}{\bar{A}}_{j}^{*}\\ 0 & -{\gamma }_{j} & {\rm{i}}{\kappa }_{j}{\bar{A}}_{i}^{*}\\ {\rm{i}}{\kappa }_{S}{\bar{A}}_{j} & {\rm{i}}{\kappa }_{S}{\bar{A}}_{i} & -{\gamma }_{S}\end{array}\right)\left(\begin{array}{c}{\scr{\Delta}} {A}_{i}\\ {\scr{\Delta}} {A}_{j}\\ {\scr{\Delta}} {A}_{S}\end{array}\right)\\ & & +\;\left(\begin{array}{ccc}0 & {\rm{i}}{\kappa }_{i}{\bar{A}}_{S} & 0\\ {\rm{i}}{\kappa }_{j}{\bar{A}}_{S} & 0 & 0\\ 0 & 0 & 0\end{array}\right)\left(\begin{array}{c}{\scr{\Delta}} {A}_{i}^{*}\\ {\scr{\Delta}} {A}_{j}^{*}\\ {\scr{\Delta}} {A}_{S}^{*}\end{array}\right)+{\rm{i}}\left(\begin{array}{c}{\gamma }_{i}{\chi }_{i}{F}_{i}\\ {\gamma }_{j}{\chi }_{j}{F}_{j}\\ {\gamma }_{S}{\chi }_{S}{F}_{S}\end{array}\right),\end{array}\end{eqnarray*}$$

where we have defined coupling parameters

$$\begin{eqnarray}&&{\kappa }_{k}={\displaystyle \frac{g{\gamma }_{k}{\chi }_{k}}{2}}={\displaystyle \frac{g}{2{m}_{k}{\omega }_{k}}};\quad k\in [i,j,S]\end{eqnarray} \tag{ A.1 }$$

for notational simplicity. The thermomechanical noise forces are assumed to be white noise correlated and obey

$$\begin{eqnarray}&&\langle {F}_{i}(t)\rangle =\langle {F}_{i}(t){F}_{j}(t\prime )\rangle =0,\end{eqnarray} \tag{ A.2 }$$

$$\begin{eqnarray}&&\langle {F}_{i}(t){F}_{j}^{*}(t+\tau )\rangle =8{\gamma }_{i}{m}_{i}{k}_{B}T{{\scr{\Delta}} }_{{ij}}{\scr{\Delta}} (\tau ).\end{eqnarray} \tag{ A.3 }$$

We decompose the complex displacements into real quadratures according to ${\scr{\Delta}} {\bf{A}}={\scr{\Delta}} \mathop{\alpha }\limits^{\rightarrow }+{\rm{i}}{\scr{\Delta}} \mathop{\beta }\limits^{\rightarrow }$ and where ${\scr{\Delta}} {\bf{A}}={({\scr{\Delta}} {A}_{i},{\scr{\Delta}} {A}_{j},{\scr{\Delta}} {A}_{S})}^{T}$. Correspondingly the noise, ${\bf{v}}={\rm{i}}{({\gamma }_{i}{\chi }_{i}{F}_{i},{\gamma }_{j}{\chi }_{j}{F}_{j},{\gamma }_{S}{\chi }_{S}{F}_{S})}^{T}\mathrm{is}$ also decomposed into real and imaginary parts, ${\bf{v}}={{\bf{v}}}_{\alpha }+{\rm{i}}{{\bf{v}}}_{\beta }$ . Expressing equation (A.1) in terms of these quantities gives

$$\begin{eqnarray}&&{\scr{\Delta}} \stackrel{\dot{}}{\mathop{\alpha }\limits^{\rightarrow }}={{\bf{M}}}_{\alpha }{\scr{\Delta}} \mathop{\alpha }\limits^{\rightarrow }+{{\bf{v}}}_{\alpha }\end{eqnarray} \tag{ A.4 }$$

$$\begin{eqnarray}&&{\scr{\Delta}} \stackrel{\dot{}}{\mathop{\beta }\limits^{\rightarrow }}={{\bf{M}}}_{\beta }{\scr{\Delta}} \mathop{\beta }\limits^{\rightarrow }+{{\bf{v}}}_{\beta }.\end{eqnarray} \tag{ A.5 }$$

For the general case, valid both above and below threshold

$$\begin{eqnarray}{{\bf{M}}}_{\alpha /\beta }={\displaystyle \frac{1}{2}}\left(\begin{array}{ccc}-{\gamma }_{i} & \mp {\kappa }_{i}| {\bar{A}}_{S}| & {\kappa }_{i}| {\bar{A}}_{j}| \\ \mp {\kappa }_{j}| {\bar{A}}_{S}| & -{\gamma }_{j} & {\kappa }_{j}| {\bar{A}}_{i}| \\ -{\kappa }_{S}| {\bar{A}}_{j}| & -{\kappa }_{S}| {\bar{A}}_{i}| & -{\gamma }_{S}\end{array}\right)\end{eqnarray} \tag{ A.6 }$$

and the elements of ${{\bf{v}}}_{\alpha ,\beta }$ satisfy $\langle {v}_{i}\rangle =0$, $\langle {v}_{k}(t){v}_{l}(t+\tau )\rangle =\frac{{\gamma }_{l}{k}_{B}T}{{m}_{l}{\omega }_{l}^{2}}{{\scr{\Delta}} }_{{kl}}{\scr{\Delta}} (\tau )$. In writing equation (A.6), we have made a choice for the pump drive phase (${\phi }_{S}=0$) and the resonator mode detection phases (${\phi }_{i,j}=0$). In general, these phases can also be chosen such that there is a coupling between the ${\scr{\Delta}} \mathop{\alpha }\limits^{\rightarrow }$ and ${\scr{\Delta}} \mathop{\beta }\limits^{\rightarrow }$ quadratures. When the pump drive phase is not fixed, but evolving in time, for instance with the pump drive being detuned, this coupling is physical and cannot be made to vanish through a suitable choice of detection phases.

The noise spectral density in the steady state is obtained by taking the expectation value after Fourier transforming and inverting equations (A.4) and (A.5), and are given by the matrix equation

$$\begin{eqnarray}&&{{\bf{S}}}_{\alpha /\beta }(\omega )={\displaystyle \frac{1}{2\pi }}{({{\bf{M}}}_{\alpha /\beta }+{\rm{i}}\omega {\bf{I}})}^{-1}{\bf{D}}{({{\bf{M}}}_{\alpha /\beta }^{T}-{\rm{i}}\omega {\bf{I}})}^{-1},\end{eqnarray} \tag{ A.7 }$$

where ${\bf{I}}$ is the identity and

$$\begin{eqnarray}{\bf{D}}=\langle {\bf{v}}{{\bf{v}}}^{T}\rangle ={k}_{B}T\left(\begin{array}{ccc}\frac{{\gamma }_{i}}{{m}_{i}{\omega }_{i}^{2}} & 0 & 0\\ 0 & \frac{{\gamma }_{j}}{{m}_{j}{\omega }_{j}^{2}} & 0\\ 0 & 0 & \frac{{\gamma }_{S}}{{M}_{S}{\omega }_{S}^{2}}\end{array}\right)\end{eqnarray} \tag{ A.8 }$$

is a matrix characterizing the diffusion due to thermal forces. The variances in steady state can be obtained from the spectrum using the Wiener-Khintchine theorem, by integrating the fluctuations over frequency, i.e.

$$\begin{eqnarray}&&{\boldsymbol{\sigma }}_{\alpha /\beta }={{\displaystyle \int }}_{-\infty }^{\infty }{{\bf{S}}}_{\alpha /\beta }(\omega ){\rm{d}}\omega .\end{eqnarray} \tag{ A.9 }$$

For finite pump detuning below threshold, the actuating force is a slowly varying function of time, i.e. ${\tilde{F}}_{S}=| {\tilde{F}}_{S}| {\rm{e}}^{-{\rm{i}}{\scr{\Delta}} t}$ where the drive frequency ${\omega }_{d}={\omega }_{S}+{\scr{\Delta}}$. This results in a pump amplitude given by ${\bar{A}}_{S}={\rm{i}}| {\bar{A}}_{S}| {\rm{e}}^{-{\rm{i}}{\scr{\Delta}} t}$.

Linearizing about the steady state amplitude, i.e ${A}_{k}={\bar{A}}_{k}+{\scr{\Delta}} {A}_{k}(t),$ where $k\in [i,j,S]$, with ${\bar{A}}_{i,j}=0$, and defining the vectors ${\scr{\Delta}} \mathop{A }\limits^{\rightarrow }={({\scr{\Delta}} {A}_{i},{\scr{\Delta}} {A}_{j})}^{T}$ and ${\scr{\Delta}} \mathop{v}\limits^{\rightarrow }={({v}_{i},{v}_{j})}^{T}$, the relevant equations of motion for the fluctuations of the resonator modes reduce to

$$\begin{eqnarray}\begin{array}{l}2\stackrel{\dot{}}{{\scr{\Delta}} \mathop{A }\limits^{\rightarrow }} & = & -\left(\begin{array}{cc}{\gamma }_{i} & 0\\ 0 & {\gamma }_{j}\ \end{array}\right){\scr{\Delta}} \mathop{A }\limits^{\rightarrow }\\ & & -\;\left(\begin{array}{cc}0 & {\kappa }_{i}| {\bar{A}}_{S}| {\rm{e}}^{-{\rm{i}}{\scr{\Delta}} t}\\ {\kappa }_{j}| {\bar{A}}_{S}| {\rm{e}}^{-{\rm{i}}{\scr{\Delta}} t} & 0\ \end{array}\right){\scr{\Delta}} {\mathop{A }\limits^{\rightarrow }}^{*}+2\mathop{v}\limits^{\rightarrow }.\end{array}\end{eqnarray} \tag{ B.1 }$$

By going to a frame rotating at $\frac{{\scr{\Delta}} }{2}$, we rewrite ${\scr{\Delta}} \mathop{A }\limits^{\rightarrow }={\scr{\Delta}} \mathop{B}\limits^{\rightarrow }{\rm{e}}^{{\displaystyle \frac{{\rm{i}}{\scr{\Delta}} t}{2}}}$ in terms of which equation (B.1) becomes

$$\begin{eqnarray}\begin{array}{l}2\stackrel{\dot{}}{{\scr{\Delta}} \mathop{B}\limits^{\rightarrow }} & = & -\left(\begin{array}{cc}{\gamma }_{i} & 0\\ 0 & {\gamma }_{j}\ \end{array}\right){\scr{\Delta}} \mathop{B}\limits^{\rightarrow }-{\rm{i}}{\scr{\Delta}} {\scr{\Delta}} \mathop{B}\limits^{\rightarrow }\\ & & -\;\left(\begin{array}{cc}0 & {\kappa }_{i}| {\bar{A}}_{S}| \\ {\kappa }_{j}| {\bar{A}}_{S}| & 0\end{array}\right){\scr{\Delta}} {\mathop{B}\limits^{\rightarrow }}^{*}+2\mathop{v}\limits^{\rightarrow }{e}^{-{\displaystyle \frac{{\scr{\Delta}} t}{2}}},\end{array}\end{eqnarray} \tag{ B.2 }$$

where ${\scr{\Delta}} \mathop{B}\limits^{\rightarrow }$ are the complex amplitudes of motion, measured at frequencies that are detuned from the individual mechanical modes ${\omega }_{i,j}$ by $\frac{{\scr{\Delta}} }{2}$. We rewrite the complex amplitudes in terms of the real quadratures ${\alpha }_{i,j},{\beta }_{i,j}$ and decompose the noise term into real and imaginary parts, i.e. ${\scr{\Delta}} \mathop{B}\limits^{\rightarrow }={\scr{\Delta}} \mathop{\alpha }\limits^{\rightarrow }+{\rm{i}}{\scr{\Delta}} \mathop{\beta }\limits^{\rightarrow }$ and $\mathop{v}\limits^{\rightarrow }={\mathop{v}\limits^{\rightarrow }}_{\alpha }+{\rm{i}}{\mathop{v}\limits^{\rightarrow }}_{\beta }$, in terms of which the equations of motion become

$$\begin{eqnarray}&&{\scr{\Delta}} \stackrel{\dot{}}{\mathop{\alpha }\limits^{\rightarrow }}={{\bf{M}}}_{\alpha }{\scr{\Delta}} \mathop{\alpha }\limits^{\rightarrow }-{\displaystyle \frac{{\scr{\Delta}} }{2}}{\scr{\Delta}} \mathop{\beta }\limits^{\rightarrow }+{\mathop{v}\limits^{\rightarrow }}_{\alpha }\end{eqnarray} \tag{ B.3 }$$

$$\begin{eqnarray}&&{\scr{\Delta}} \stackrel{\dot{}}{\mathop{\beta }\limits^{\rightarrow }}={{\bf{M}}}_{\beta }{\scr{\Delta}} \mathop{\beta }\limits^{\rightarrow }+{\displaystyle \frac{{\scr{\Delta}} }{2}}{\scr{\Delta}} \mathop{\alpha }\limits^{\rightarrow }+{\mathop{v}\limits^{\rightarrow }}_{\beta },\end{eqnarray} \tag{ B.4 }$$

where

$$\begin{eqnarray}{{\bf{M}}}_{\alpha /\beta }={\displaystyle \frac{1}{2}}\left(\begin{array}{cc}-{\gamma }_{i} & \mp {\kappa }_{i}| {\bar{A}}_{S}| \\ \mp {\kappa }_{j}| {\bar{A}}_{S}| & -{\gamma }_{j}\end{array}\right)\end{eqnarray} \tag{ B.5 }$$

and the elements of ${\mathop{v}\limits^{\rightarrow }}_{\alpha ,\beta }$ satisfy $\langle {v}_{k,\eta }\rangle =0,\eta \in [\alpha ,\beta ];k\in [i,j]$ and $\langle {v}_{k,\eta }(t){v}_{l,{\eta }^{\prime }}(t+\tau )\rangle ={\displaystyle \frac{{\gamma }_{l}{k}_{B}T}{{m}_{l}{\omega }_{l}^{2}}}{{\scr{\Delta}} }_{{kl}}{{\scr{\Delta}} }_{\eta ,{\eta }^{\prime }}{\scr{\Delta}} (\tau )$.

The coupling between the ${\scr{\Delta}} \mathop{\alpha }\limits^{\rightarrow }$ and ${\scr{\Delta}} \mathop{\beta }\limits^{\rightarrow }$ quadratures of the individual oscillators resulting from the detuned drive is apparent in the above equations. The steady state correlations between these quadratures can be obtained by forming the following four-dimensional vectors; ${\bf{Z}}={({\scr{\Delta}} {\alpha }_{i},{\scr{\Delta}} {\alpha }_{j},{\scr{\Delta}} {\beta }_{i},{\scr{\Delta}} {\beta }_{j})}^{T}={({\scr{\Delta}} \mathop{\alpha }\limits^{\rightarrow },{\scr{\Delta}} \mathop{\beta }\limits^{\rightarrow })}^{T}$ and ${\bf{v}}={({\mathop{v}\limits^{\rightarrow }}_{\alpha },{\mathop{v}\limits^{\rightarrow }}_{\beta })}^{T}$, in terms of which the equations of motion become

$$\begin{eqnarray}&&\stackrel{\dot{}}{{\bf{Z}}}={\bf{M}}{\bf{Z}}+{\bf{v}}\end{eqnarray} \tag{ B.6 }$$

$$\begin{eqnarray}{\bf{M}}=\left(\begin{array}{cc}{{\bf{M}}}_{\alpha } & -\frac{{\scr{\Delta}} }{2}{\bf{I}}\\ \frac{{\scr{\Delta}} }{2}{\bf{I}} & {{\bf{M}}}_{\beta }\end{array}\right)\end{eqnarray} \tag{ B.7 }$$

and ${\bf{I}}$ is the $2\times 2$ identity matrix.

The noise spectral densities are obtained by solving equation (B.6) in Fourier space, as before. The spectrum in the steady state is

$$\begin{eqnarray}\begin{array}{l}{\bf{S}}(\omega ) & = & \left\langle{\bf{Z}}(\omega ){\bf{Z}}{(\omega )}^{\dagger }\right\rangle\\ & = & {\displaystyle \frac{1}{2\pi }}{({\bf{M}}+{\rm{i}}\omega {\bf{I}})}^{-1}{\bf{D}}{({{\bf{M}}}^{T}-{\rm{i}}\omega {\bf{I}})}^{-1},\end{array}\end{eqnarray} \tag{ B.8 }$$

where ${\bf{I}}$ is the identity and

$$\begin{eqnarray}{\bf{D}}=\langle {\bf{v}}{{\bf{v}}}^{T}\rangle ={k}_{B}T\left(\begin{array}{cccc}\frac{{\gamma }_{i}}{{m}_{i}{\omega }_{i}^{2}} & 0 & 0 & 0\\ 0 & \frac{{\gamma }_{j}}{{m}_{j}{\omega }_{j}^{2}} & 0 & 0\\ 0 & 0 & \frac{{\gamma }_{i}}{{m}_{i}{\omega }_{i}^{2}} & 0\\ 0 & 0 & 0 & \frac{{\gamma }_{j}}{{m}_{j}{\omega }_{j}^{2}}\end{array}\right).\end{eqnarray} \tag{ B.9 }$$

We construct composite quadratures ${x}_{\pm }=({\alpha }_{i}\pm {\alpha }_{j})/\sqrt{2}$, ${y}_{\pm }=({\beta }_{i}\pm {\beta }_{j})/\sqrt{2}$ as before and represent the fluctuations in these quadratures by the column matrix, ${{\bf{Z}}}_{c}={({\scr{\Delta}} {x}_{+},{\scr{\Delta}} {x}_{-},{\scr{\Delta}} {y}_{+},{\scr{\Delta}} {y}_{-})}^{T}$, which is related to ${\bf{Z}}$ by

$$\begin{eqnarray}{{\bf{Z}}}_{c}={\bf{R}}{\bf{Z}};\quad {\bf{R}}={\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{cccc}1 & 1 & 0 & 0\\ 1 & -1 & 0 & 0\\ 0 & 0 & 1 & 1\\ 0 & 0 & 1 & -1\end{array}\right).\end{eqnarray} \tag{ B.10 }$$

The correlations of the composite quadratures are therefore given by

$$\begin{eqnarray}&&{{\bf{S}}}_{c}(\omega )=\left\langle{{\bf{Z}}}_{c}(\omega ){{\bf{Z}}}_{c}{(\omega )}^{\dagger }\right\rangle={\bf{R}}{\bf{S}}{{\bf{R}}}^{T}.\end{eqnarray} \tag{ B.11 }$$

We consider the case where the frequencies of the resonator modes are identical and the losses are symmetric (${{\scr{\Delta}} }_{\gamma }={{\scr{\Delta}} }_{\omega }=0$). For this case, the diffusion matrix ${\bf{D}}=\frac{{k}_{B}T\gamma }{m{\omega }^{2}}{{\scr{\Delta}} }_{{ij}}$.

The correlations between the composite quadratures in this case are given by

$$\begin{eqnarray}{{\bf{S}}}_{c}(\omega )=\left(\begin{array}{cccc}{S}_{{x}_{+},{x}_{+}} & 0 & {S}_{{x}_{+},{y}_{+}} & 0\\ 0 & {S}_{{x}_{-},{x}_{-}} & 0 & {S}_{{x}_{-},{y}_{-}}\\ {S}_{{y}_{+},{x}_{+}} & 0 & {S}_{{y}_{+},{y}_{+}} & 0\\ 0 & {S}_{{y}_{-},{x}_{-}} & 0 & {S}_{{y}_{-},{y}_{-}}\end{array}\right),\end{eqnarray} \tag{ B.12 }$$

where the non zero correlations are as indicated above. There are no correlations between (${x}_{+}$, ${x}_{-}$) and (${y}_{+}$, ${y}_{-}$), given our choice of detection phases and the fact that we consider the case ${{\scr{\Delta}} }_{\gamma }={{\scr{\Delta}} }_{\omega }=0$. The correlations in equation (B.12) evaluate to

$$\begin{eqnarray}\begin{array}{l}{S}_{{x}_{\pm },{x}_{\pm }}(\omega ) & = & {S}_{{y}_{\mp },{y}_{\mp }}(\omega )\\ & = & {\displaystyle \frac{{k}_{B}T\gamma }{m{\omega }^{2}}}\;{\displaystyle \frac{2\left({{\scr{\Delta}} }^{2}+{\gamma }^{2}{(1\mp \mu )}^{2}+4{\omega }^{2}\right)}{\pi \left(4{\omega }^{2}+{\lambda }_{+}^{2}\right)\left(4{\omega }^{2}+{\lambda }_{-}^{2}\right)}}\end{array}\end{eqnarray} \tag{ B.13 }$$

$$\begin{eqnarray}\begin{array}{l}{S}_{{x}_{+},{y}_{+}}(\omega ) & = & {S}_{{y}_{+},{x}_{+}}(\omega )\\ & = & {\displaystyle \frac{{k}_{B}T\gamma }{m{\omega }^{2}}}\;{\displaystyle \frac{4{\scr{\Delta}} (\gamma \mu +2{\rm{i}}\omega )}{\pi \left(4{\omega }^{2}+{\lambda }_{+}^{2}\right)\left(4{\omega }^{2}+{\lambda }_{-}^{2}\right)}}\\ & = & {S}_{{x}_{-},{y}_{-}}(-\omega )={S}_{{y}_{-},{x}_{-}}(-\omega ),\end{array}\end{eqnarray} \tag{ B.14 }$$

where ${\lambda }_{\pm }^{2}={\gamma }^{2}\left(1+{\mu }^{2}\right)-{{\scr{\Delta}} }^{2}\pm 2\gamma \sqrt{{\gamma }^{2}{\mu }^{2}-{{\scr{\Delta}} }^{2}}$. The variances in the steady state obtained using the Weiner-Khintchine theorem are

$$\begin{eqnarray}\begin{array}{l}{\sigma }_{{x}_{\pm },{x}_{\pm }} & = & \left({\displaystyle \frac{{k}_{B}T\gamma }{m{\omega }^{2}}}\right)\left[{\displaystyle \frac{\left({{\scr{\Delta}} }^{2}+{\gamma }^{2}{(1\mp \mu )}^{2}-{\lambda }_{-}^{2}\right)}{{\lambda }_{+}{\lambda }_{-}\left({\lambda }_{+}+{\lambda }_{-}\right)}}+{\displaystyle \frac{1}{{\lambda }_{+}}}\right]\\ & = & {\sigma }_{{y}_{\mp },{y}_{\mp }}.\end{array}\end{eqnarray} \tag{ B.15 }$$

Unlike the case of zero detuning, we obtain non-zero steady state correlations between the x and y quadratures

$$\begin{eqnarray}&&{\sigma }_{{x}_{\pm },{y}_{\pm }}=\left({\displaystyle \frac{{k}_{B}T\gamma }{m{\omega }^{2}}}\right){\displaystyle \frac{2{\scr{\Delta}} \gamma \mu }{{\lambda }_{+}{\lambda }_{-}\left({\lambda }_{+}+{\lambda }_{-}\right)}}={\sigma }_{{y}_{\pm },{x}_{\pm }}.\end{eqnarray} \tag{ B.16 }$$

The correlations between the ${x}_{\pm }$ and ${y}_{\pm }$ are now proportional to the drive detuning.