Mechanical entanglement detection in an optomechanical system

TL;DR

This paper introduces a novel optomechanical setup with multiple probing tones to generate and detect entanglement between two mechanical resonators by evaluating their collective quadratures.

Contribution

It presents a new scheme combining two driving pumps and four probing tones for entanglement generation and detection in mechanical resonators.

Findings

Enables generation of two-mode mechanical squeezing.

Allows separate evaluation of collective mechanical quadratures.

Facilitates inference of entanglement between resonators.

Abstract

We propose here a setup to generate and evaluate the entanglement between two mechanical resonators in a cavity optomechanical setting. As in previous proposals, our scheme includes two driving pumps allowing for the generation of two-mode mechanical squeezing. In addition, we include here four additional probing tones, which allow for the separate evaluation of the collective mechanical quadratures required to estimate the Duan quantity, thus allowing us to infer whether the mechanical resonators are entangled.

Figures

Equations

$$⟨\Delta X_{\Sigma }^2⟩+⟨\Delta Y_{\Delta }^2⟩\le 1$$

$$H={\omega }_{\mathrm{c}}{a}^{†}a+$$

$$+$$

$$H=$$

$$+$$

$${\chi }_{\mathrm{c}}^{-1}\left(\omega +\delta \right)a_{\omega }=$$

$$-i{\mathcal{G}}_{\Delta 1}\left[{\beta }_{1\omega -\delta }+{\beta }_{2\omega -\delta }\right]$$

$$-i{\mathcal{G}}_{\Delta 2}\left[{\beta }_{1\omega +\delta }+{\beta }_{2\omega +\delta }\right]+\sqrt{\kappa }{a}_{\omega }^{\mathrm{in}}$$

$${\chi }_{\mathrm{m}}^{-1}\left(\omega -\delta \right){\beta }_{\omega ,1}=$$

$$-i{\mathcal{G}}_{\Delta 2}^{\ast }{a}_{\omega +\delta }+\sqrt{\gamma }{\beta }_{1\omega }^{\mathrm{in}}$$

$${\chi }_{\mathrm{m}}^{-1}\left(\omega +\delta \right){\beta }_{2\omega }=$$

$$-i{\mathcal{G}}_{\Delta 2}^{\ast }{a}_{\omega +\delta }+\sqrt{\gamma }{\beta }_{2\omega }^{\mathrm{in}}$$

$$\sqrt{\kappa }{a}_{\omega }^{\mathrm{in}}=\sqrt{{\kappa }_{\mathrm{i}}}{a}_{\mathrm{i}\omega }^{\mathrm{in}}+\sqrt{{\kappa }_{\mathrm{e}}}{a}_{\mathrm{e}\omega }^{\mathrm{in}}.$$

$${\beta }_{1\omega }=u{b}_{1\omega }+v{b}_{2-\omega }^{†},$$

$${\beta }_{2\omega }=u{b}_{2\omega }+v{b}_{1-\omega }^{†},$$

$$S(r)=\exp \left(r{b}_{\omega ,1}{b}_{\omega ,2}-h.c.\right)$$

$$S_{\omega }^{\Sigma \theta }=\frac{1}{2}⟨\left\{⟨X_{-\omega }^{\Sigma \theta },X_{\omega }^{\Sigma \theta }⟩\right\}⟩$$

$$S_{\omega }^{\Delta \theta }=\frac{1}{2}⟨\left\{⟨X_{-\omega }^{\Delta \theta },X_{\omega }^{\Delta \theta }⟩\right\}⟩$$

$$X_{\omega }^{\Sigma }\doteq X_{1\omega }+X_{2\omega }$$

$$X_{\omega }^{\Delta }\doteq X_{1\omega }-X_{2\omega }$$

$$X_{\omega }^{\mathrm{out}}$$

$$\sqrt{\kappa }$$

$$+(u-v)$$

$${\bar{X}}_{\omega }^{\Sigma }\doteq {\bar{X}}_{1\omega }+{\bar{X}}_{2\omega }$$

$${\bar{X}}_{\omega }^{\Delta }\doteq {\bar{X}}_{1\omega }-{\bar{X}}_{2\omega }$$

$${\bar{X}}_{1\omega }$$

$${\bar{X}}_{2\omega }$$

$${\bar{X}}_{1t}=b_{\mathrm{t}}{e}^{i\delta t}+b_{\mathrm{t}}^{†}{e}^{-i\delta t}$$

$$⟨\Delta {\bar{X}}_{\Sigma }^2⟩+⟨\Delta {\bar{Y}}_{\Delta }^2⟩=⟨\Delta X_{\Sigma }^2⟩+⟨\Delta Y_{\Delta }^2⟩,$$

$$S_{\omega }^{\mathrm{out}\theta }={\mathcal{G}}_D^2{\left|{\chi }_{\mathrm{c}}(\delta )\right|}^2$$

$$\{{\left(u+v\right)}^2\left[{\left(\cos \theta \cos \phi \right)}^2{\bar{S}}_{\omega }^{\Sigma 0}+{\left(\cos \theta \sin \phi \right)}^2{\bar{S}}_{\omega }^{\Delta \pi /2}\right]$$

$$+{\left(u-v\right)}^2\left[{\left(\sin \theta \cos \phi \right)}^2{\bar{S}}_{\omega }^{\Sigma \pi /2}+{\left(\sin \theta \sin \phi \right)}^2{\bar{S}}_{\omega }^{\Delta 0}\right]\}$$

$$+B^{\mathrm{in}}+C_{\omega }^{\mathrm{in}}\cos \left[2\varphi \right]$$

$${S_{\omega }^{\mathrm{out}0}|}_{\phi =0}$$

$${S_{\omega }^{\mathrm{out}0}|}_{\phi =\pi /2}$$

$${S_{\omega }^{\mathrm{out}\pi /2}|}_{\phi =0}$$

$${S_{\omega }^{\mathrm{out}\pi /2}|}_{\phi =\pi /2}$$

$$⟨\Delta {\bar{X}}_{\Sigma }^2⟩+⟨\Delta {\bar{Y}}_{\Delta }^2⟩=$$

$$\left[\gamma \kappa \left(n_{\mathrm{m}}+\frac{1}{2}\right)+4\left(G_{-}^2+G_{+}^2\right)\left(n_{\mathrm{c}}+\frac{1}{2}\right)\right]/\left[2(G_{-}^2-G_{+}^2)\right]$$

$$H=$$

Full text

Since the early years of quantum mechanics, it was realised that some of the consequences borne from its fundamental principles are in stark contradiction with an intuitive picture of reality deriving from our daily experience of the world. Arguably, one of its most unsettling aspects–yet potentially useful in the manipulation of information at the quantum level– is represented by entanglement, which is a form of correlation inherent to quantum mechanical systems. One of the prototypical examples of entangled systems was discussed as early as 1935 by Einstein Podolsky and Rosen Einstein et al. (1935), in the attempt to prove the incompleteness of quantum mechanics. In their paper the authors discuss a gedankenexperiment in which they show how the measurement of position or momentum on a quantum system can affect the state of a second system causally disconnected from the first one. While in the decades following the paper by Einstein Podolsky and Rosen, the reality of entanglement has been demonstrated in the context of quantum optics (see e.g. Aspect et al. (1982); Reid and Walls (1986); Pan et al. (2003)), the realisation of the experiment involving the position and momentum of a (macroscopic) material system, much along the lines of the original proposal discussed in Einstein et al. (1935), has not been carried out.

In the recent years the progress in the physics of optomechanical systems Milburn and Woolley (2012); Aspelmeyer et al. (2014) has opened the prospect of entangling the mechanical degrees of freedom of two macroscopic mechanical oscillators Mancini et al. (2002); Pinard et al. (2005); Pirandola et al. (2006); Schmidt et al. (2012); Abdi et al. (2012); Tan et al. (2013); Wang and Clerk (2013); Woolley and Clerk (2014); Szorkovszky et al. (2014), following the experimental realisation entanglement between a mechanical oscillator with a microwave field Palomaki et al. (2013) and the backaction-evading measurement of collective mechanical modes Ockeloen-Korppi et al. (2016). In particular in Wang and Clerk (2013); Woolley and Clerk (2014); Li et al. (2015) it has been proposed to use a common optical cavity with an appropriate optical drive in order to entangle the degrees of freedom of two mechanical resonators. While the realisation of such a setup is within reach of the current experimental capabilities Massel et al. (2012); Ockeloen-Korppi et al. (2016), a reliable characterisation of the entanglement properties between mechanical resonators, considering the limitations imposed by the current experimental parameters range, is lacking. In particular, in the theoretical proposals based on a one-cavity/two-mechanical resonators concept, the most severe limitation is represented by the impossibility, in the current experimental settings, of separately addressing the two mechanical modes. The separate addressability of the mechanical modes, which, in principle, could be realised by engineering mechanical resonators of sufficiently different frequency, would open up the prospective of state tomography for the two mechanical resonators Vanner et al. (2015) and thus provide a different route to the characterisation of the entanglement properties of the two mechanical resonators.

In this article, we propose a detection scheme for currently available one-cavity/two-resonators setups related to the concept of backaction-evading measurement Caves et al. (1980); Braginsky et al. (1980); Clerk et al. (2008, 2010); Hertzberg et al. (2010); Suh et al. (2014); Ockeloen-Korppi et al. (2016), in a system constituted by one resonant cavity coupled to two mechanical resonators by radiation pressure force (see Fig. 1). The scheme proposed here, based on a novel pump/probe setup, allows us to verify the entanglement between the two mechanical resonators.

As previously noted in the literature Woolley and Clerk (2014), it is possible to induce two-mode squeezing Walls and Milburn (2008) on the mechanical degrees of freedom by suitably driving the cavity with two coherent tones. While this two-mode squeezing represents the key element in the definition of the entanglement between the two mechanical operators, the direct observation of entanglement between the mechanical resonators has proven elusive. Here we discuss how, in order to gain access to the different quadratures of the collective modes needed to determine the entanglement between mechanical resonators, 4 extra probes are required. More specifically, the setup proposed here allows us to infer the value of the Duan quantity Duan et al. (2000) for the collective quadratures associated with the dynamics of the two mechanical resonators from the noise spectrum of the output field.

We evaluate here the Duan quantity and the ensuing Duan bound Duan et al. (2000), which represent a possible inseparability criterion for a bipartite system and have been previously considered in the context of entanglement between mechanica resonators in cavity optomechanics Pinard et al. (2005); Woolley and Clerk (2014); Li et al. (2015).

The Duan bound is expressed in terms of an upper bound to the total variance of EPR-like observables. For instance, for two mechanical resonators, the Duan criterion could be expressed as

[TABLE]

where $X_{\Sigma}=X_{1}+X_{2}$ and $Y_{\Delta}=Y_{1}-Y_{2}$, with $X_{1,2}$ and $Y_{1,2}$ representing two orthogonal quadrature operators( e.g. the position and momentum operators) associated with the dynamics of subsystem $1$ and $2$ respectively.

Setup and equations of motion–

The setup considered here consists of two mechanical resonators dispersively coupled to a single optical cavity. The general Hamiltonian of the (isolated) system can be written as

[TABLE]

where $a$, $b_{1}$ and $b_{2}$, along with their hermitian conjugates, represent the field operators associated with the cavity and the two mechanical resonators, with resonant frequencies $\omega_{\mathrm{c}}$, $\omega_{1}$ and $\omega_{2}$ respectively. Furthermore, we have assumed that each mechanical resonator is coupled through a radiation-pressure coupling term of strength $g_{0}$. In our scheme the system is strongly driven at frequencies $\omega_{\mathrm{c}}+\omega_{1}$ with a coherent pumping tone of amplitude $\alpha_{+}$ and $\omega_{\mathrm{c}}-\omega_{2}$ with amplitude $\alpha_{-}$ ($G_{\pm}=g_{0}\alpha_{\pm}$). In addition to the driving at $\omega_{\mathrm{c}}\pm\omega_{1,2}$, the cavity driving scheme proposed here comprises four extra detection tones at frequencies $\omega_{\mathrm{c}}+\omega_{1}\pm\delta$ and $\omega_{\mathrm{c}}-\omega_{2}\pm\delta$, with $\delta=\left(\omega_{1}-\omega_{2}\right)/2$ (see Fig. 1).

In the appropriate frame ($\omega_{\mathrm{c}}+\delta$ and $\omega_{\Sigma}=(\omega_{1}+\omega_{2})/2$ for cavity and mechanics respectively), after linearisation around the pumping tones, and neglecting terms rotating at frequencies $\omega\simeq 2\omega_{\Sigma}$ since we assume $\omega_{\Sigma}\gg\kappa$ (rotating-wave approximation), the Hamiltonian can be written as

[TABLE]

with $G_{p,q\pm}=g_{0}\alpha_{p,q\pm}$.

Choosing $G_{+}/G_{-}=G_{p+}/G_{p-}=G_{q+}/G_{q-}$, we can write the quantum Langevin equations (QLEs) associated with the linearized system Hamiltonian (2), in the Fourier domain, as (see A)

[TABLE]

where $\mathcal{G}=\sqrt{\left|G_{-}\right|^{2}-\left|G_{+}\right|^{2}}$, $\mathcal{G}_{\Delta 1}=\left[G_{-}G_{p-}^{*}-G_{+}G_{p+}^{*}\right]/\mathcal{G}$, $\mathcal{G}_{\Delta 2}=\left[G_{-}G_{q-}^{*}-G_{+}G_{q+}^{*}\right]/\mathcal{G}$ and $\chi_{\mathrm{x}}(\omega)^{-1}=\gamma_{\mathrm{x}}/2-i\omega$ ($\mathrm{x}=1,2$, $\gamma_{1}=\gamma$,$\gamma_{2}=\kappa$). The input field introduced in (3) include contributions both from internal and external noise (see Fig. 1)

[TABLE]

Most importantly, the mechanical operators

[TABLE]

where $u=G_{+}/\sqrt{G_{+}^{2}-G_{-}^{2}}$, $v=G_{-}/\sqrt{G_{+}^{2}-G_{-}^{2}}$, can be regarded as being generated by the action of the two-mode squeezing operator

[TABLE]

with $r=\operatorname{arccosh}u$. $S(r)$ can be shown to give rise to quantum correlations between the mechanical modes, analogously to the situation encountered in the context of quantum optics, e.g. in the case of non-degenerate parametric amplification Walls and Milburn (2008) and thus represents the key ingredient for entanglement generation in the present setup Woolley and Clerk (2014).

Our strategy in the solution of the problem in the presence of both pumping and probing tones, consists now in assuming that $\left|\mathcal{G}\right|\gg\left|\mathcal{G}_{\Delta 1}\right|,\left|\mathcal{G}_{\Delta 2}\right|$, and $\gamma\ll\delta$, treating the probing tones as a perturbation with respect to the driving tones (see B). These conditions allow us to solve for the dynamics of the mechanical resonators as if it was determined by the pumping tone only and independently for each mode.

In Fig. 2, we have depicted the mechanical noise spectrum, with the following definitions

[TABLE]

and

[TABLE]

The spectra in Fig. 2 are obtained in the presence of thermal noise both for the cavity ( $\braket{a^{\mathrm{in}}_{\mathrm{I}\,\omega}{a^{in}_{\mathrm{I}\,\omega^{\prime}}}^{\dagger}}=\left(n_{\mathrm{i}}+1\right)\delta_{\omega,\omega^{\prime}}$, $\braket{a^{\mathrm{in}}_{\mathrm{E}\,\omega}{a^{in}_{\mathrm{E}\,\omega^{\prime}}}^{\dagger}}=\left(n_{\mathrm{e}}+1\right)\delta_{\omega,\omega^{\prime}}$), and the mechanical bath ($\braket{b^{\mathrm{in}}_{\mathrm{x}\,\omega}{b^{in}_{\mathrm{x}\,\omega^{\prime}}}^{\dagger}}=\left(n_{\mathrm{m}}+1\right)\delta_{\omega,\omega^{\prime}}$, $\mathrm{x}=1,2$).

The solution of the equations of motion for the mechanical degrees of freedom determines the dynamics of the cavity field, giving rise to the appearance five peaks in the cavity (and output) spectrum. Due to the small value of the mechanical linewidth, it is possible to consider separately each peak induced in the cavity field by the mechanical resonators. In our analysis, we are interested in particular in the peak at $\omega=0$ ($\omega=\omega_{\mathrm{c}}+\delta$ in the original frame), which comprises contributions from the dynamics of both mechanical resonators and, as we will show, contains all the information needed to evaluate the Duan bound. Since for $\omega\simeq 0$ the mechanical contributions to the cavity field are predominantly provided by $\mathcal{G}_{\Delta 1}\beta_{\omega-\delta,2}$ and $\mathcal{G}_{\Delta 2}\beta_{\omega+\delta,1}$, due to the resonance condition in the mechanical equations of motion: for $\omega\simeq 0$, $\beta_{\omega,1}$ is resonant at $\omega+\delta$ and $\beta_{\omega,2}$ at $\omega-\delta$.

Output spectrum–

If homodyne detection is performed on the fluctuations, from Eq. (8), expressing $\beta_{1}$ and $\beta_{2}$ in terms of $b_{1}$ and $b_{2}$, it is possible to monitor the dynamics of the collective mechanical modes through the measurement of the quadratures of the output field, namely (for $\omega\simeq 0$, $\delta\ll\kappa$)

[TABLE]

where the dynamics of the collective mechanical modes is encoded in the frequency-shifted quadrature operators $\bar{X}^{\Sigma}_{\omega}$, $\bar{Y}^{\Sigma}_{\omega}$ and $\bar{X}^{\Delta}_{\omega}$, $\bar{Y}^{\Delta}_{\omega}$ defined by

[TABLE]

with

[TABLE]

with analogous definitions holding for the quadratures $\bar{Y}_{1\,\omega}$ and $\bar{Y}_{2\,\omega}$. Note that the frequency-shifted quadrature operators defined above, do not directly correspond to the usual quadrature operators. In the time domain, the operators defined in (9) acquire a nontrivial time dependence, for instance we have that

[TABLE]

and thus correspond to the mechanical quadratures for $\delta=0$ only. While one cannot directly relate the frequency-shifted quadrature operators to the regular ones, for each pair of orthogonal quadratures, the uncertainties associated with the collective quadratures of the mechanical motion, fulfil the following relation

[TABLE]

which, crucially, allows us to establish the link between the output spectrum and the spectrum of the collective mechanical quadratures. Therefore, the knowledge of one pair of orthogonal frequency-shifted mechanical quadratures allows one to deduce the value of the corresponding pair of regular quadratures and, consequently, to infer the value of the Duan quantity from the spectrum of the output field.

In the derivation of Eq. (7), we have chosen the average phase of the detection tones as the phase reference with respect to which both the phases of the probing tones $\mathcal{G}_{\Delta_{1}}$ and $\mathcal{G}_{\Delta_{2}}$ ($\pm\varphi$) as well as the homodyne detection phase $\theta$ are referred. From Eq. (7), it is possible to relate the output spectrum $S^{\mathrm{out}\,\theta}_{\omega}$ to the spectrum of the frequency-shifted mechanical quadratures as

[TABLE]

(see the C for the derivation of the output noise spectrum $S^{\mathrm{out}\,\theta}_{\omega}$, and the definitions of $B^{\mathrm{in}}$ and $C_{\omega}^{\mathrm{in}}$).

Eq. (8) expresses the possibility to access the collective mechanical noise spectra $\bar{S}_{\omega}^{\Sigma\,0}$, $\bar{S}_{\omega}^{\Sigma\,\pi/2}$, $\bar{S}_{\omega}^{\Delta\,0}$, $\bar{S}_{\omega}^{\Delta\,\pi/2}$, by changing the relative phase of the detection tones $\varphi$ and the phase of the homodyne detector $\theta$.

The measurement strategy leading to the determination of the Duan quantity consists in the measurement of the output spectrum for four different values of $\left(\theta,\varphi\right)=\left(0,0\right)$, $\left(0,\pi/2\right)$, $\left(\pi/2,0\right)$, $\left(\pi/2,\pi/2\right)$, yielding

[TABLE]

Having in mind the relation established by Eq. (10) The Duan quantity can then be determined as the sum of the appropriate output quadrature spectra as determined in Eqs. (13). The output spectrum therefore provides a measurement of the collective mechanical quadratures induced by the pumping tones $\alpha_{\pm}$, disregarding higher-order effects induced by the detection tones. As detailed in B, for $\omega\simeq 0$ the correction to the dynamics of the mechanical resonators due to the detection tones is of the order $\left(\mathcal{G}_{D}/\mathcal{G}\right)^{2}$, and thus can be safely neglected for a suitable choice of readout tones.

In the limit $\kappa\gg\delta,\mathcal{G}\gg\gamma$, the quantity given in Eq. (10) can be written as

[TABLE]

(see C for the full exact expression), which corresponds to the approximate formula given in Woolley and Clerk (2014) for the Duan quantity in the adiabatic limit.

The detection scheme proposed here represents a somewhat idealised setup. A first apporximation in our approach is represented by the rotating-wave approximation. Considering that that discarded terms in performing the rotating wave approximation are of the order on $\sim\left(\kappa/\omega_{\mathrm{m}}\right)^{2}$, reaching a range of experimental parameters in which this approximation does not seem to represent an insurmountable experimental challenge (see e.g. Pirkkalainen et al. (2015), where $\kappa/\omega_{\mathrm{m}}\simeq 0.05$).

The most relevant source of non-ideality is probably represented by the presence of parametric modulation terms to the dynamics of the cavity. This effect was observed in the context of the optomechanical generation and detection of mechanical squeezing Wollman et al. (2015); Pirkkalainen et al. (2015) and is likely to affect the detection of mechanical entanglement as well.

Among other non-idealities likely encountered in the specific experimental realisation of the detection scheme discussed here, it is worth mentioning that the two mechanical resonators are inevitably not identical and the thermal baths to which they are coupled are not necessarily at the same temperature. The most severe constraint posed by the difference between the two mechanical resonators seems to be posed by the different value of the bare optomechanical coupling $g_{0}$, which however can be compensated by a suitably engineered pumping scheme (see e.g. Woolley and Clerk (2014)).

We have proposed here a potential scheme for the detection of entanglement between mechanical resonators coupled to a common optical cavity, establishing a straightforward setup for the quantification of such entanglement in terms of the Duan quantity. Its feasibility is within the framework of current experimental capabilities.

The author would like to thank Mika Sillanpää and Tero Heikkilä for useful discussions. This work was supported by the Academy of Finland (Contract No. 27545).

Appendix A Derivation of the equations of motion

In the appropriate frame ($\omega_{\mathrm{c}}+\delta$ and $\omega_{\Sigma}=(\omega_{1}+\omega_{2})/2$ for cavity and mechanics respectively), and invoking the rotating-wave approximation, the Hamiltonian given in Eq. (2) of the main text can be written as

[TABLE]

with $\delta=\left(\omega_{1}-\omega_{2}\right)/2$. From the Hamiltonian (15), it is possible to write down the quantum Langevin equations for cavity and mechanical degrees of freedom as

[TABLE]

where we have defined $b_{\Sigma}=\left(b_{1}+b_{2}\right)/\sqrt{2}$, $b_{\Delta}=\left(b_{1}-b_{2}\right)/\sqrt{2}$. If we now write the EOM for the linear combination $\left(G_{-}b_{\Sigma,\Delta}+G_{+}b_{\Sigma,\Delta}^{\dagger}\right)$, from Eq. (16) we can write

[TABLE]

and, analogously,

[TABLE]

The cancellations on the second and third line of Eq. (17), needed to recast the problem in terms of Bogolyubov modes, occur only if $G_{p+}/G_{p-}=G_{q+}/G_{q\,-}=G_{+}/G_{-}$. In this case,

[TABLE]

where we have adopted the definitions

[TABLE]

In order to ensure that the mechanical operators on the rhs of Eq. (16) can be all expressed in terms of the same Bogolyubov operator, we have to verify under what conditions the following three (in principle different) transformations

[TABLE]

coincide, having set

[TABLE]

It turns out that the condition $G_{p+}/G_{p-}=G_{q+}/G_{q\,-}=G_{+}/G_{-}$ is a sufficient condition to determine $u=u^{\prime}=u^{\prime\prime}$ and $v=v^{\prime}=v^{\prime\prime}$. In our analysis we have set the reference phase to be given by the pump tones $G_{+}$ and $G_{-}$. Focusing on $u^{\prime}$, we have

[TABLE]

and analogously for $v^{\prime}$, $u^{\prime\prime}$ and $v^{\prime\prime}$. Defining

[TABLE]

we can write the quantum Langevin equations of motion induced by the Hamiltonian (15) in the Fourier domain as

[TABLE]

which correspond to Eq. (3) of the main text.

Appendix B Equations of motion: perturbative solution

In the following we derive the solution of the equations of motion around $\omega\simeq 0$, treating the probing tones $\mathcal{G}_{\Delta 1}$ and $\mathcal{G}_{\Delta 2}$ perturbatively. The solution that we obtain for the output field is then used to determine the output noise spectrum, from which, in turn, it is possible to deduce the mechanical modes dynamics. The solution of the equations of motion (25), considering that $\mathcal{G}\ll\left|\mathcal{G}_{\Delta 1}\right|=\left|\mathcal{G}_{\Delta 2}\right|$, can be obtained with the aid of perturbation theory. Defining

[TABLE]

we can formally write the solution for $a_{\omega}$, $\beta_{1\,\omega}$, $\beta_{2\,\omega}$ as

[TABLE]

Substituting the perturbative expression given by Eq. (26), Eqs. (25) can be solved order-by-order in $\lambda$,

[TABLE]

The [math]-th order term of Eq. (27) is represented by the equations governing the 2-pump driving scheme, in the absence of detection tones which can be readily solved to give

[TABLE]

where

[TABLE]

The assumption $\gamma\ll\delta$ allows us to recognise that $\beta^{0}_{1\,\omega}$ and $\beta^{0}_{2\,\omega}$ are peaked around $\omega\simeq\delta$ and $\omega\simeq-\delta$ respectively, while $a_{\omega}^{0}$ exhibits a double-peak structure for $\omega\simeq\pm\delta$.

The solution given by (28), allows us to write the first-order approximation of the system equations of motion ($n=1$ in Eq. (27)). The considerations concerning the peak structure of the [math]-th order equations - essentially because $\chi_{\mathrm{m}}\left(\omega\right)\sim\delta_{\omega}$ –, allow us to write the first order approximation for the cavity field and the mechanics as

$\omega\simeq 0$

[TABLE]

The value of the first-order correction to the cavity field, can be expressed in terms of the input field

[TABLE]

On the other hand, Eq. (LABEL:eq:4x) can be shown to encode the relevant information about the mechanical quadratures

[TABLE]

where $\Phi=\phi_{1}+\phi_{2}$ and $\varphi=\phi_{1}-\phi_{2}$.

From Eq. (1) it is possible to notice that $a_{\omega}^{1}$ does not depend on the value of the input fields at $\omega=0$, implying that the backaction term, to this order, will not give rise to any interference contribution with the zeroth-order cavity field. 2. 2.

$\omega\simeq\delta$

[TABLE] 3. 3.

$\omega\simeq-\delta$

[TABLE]

The same strategy can be applied to determine the second-order approximation to the cavity field around $\omega\simeq 0$ ($n=2$ in Eq. (27)), giving

[TABLE]

and, in terms of [math]-th order approximation,

[TABLE]

As previously discussed, due to the structure of the mechanical response, the [math]-th order cavity response will be peaked around $w\simeq\pm\delta$, implying that the terms appearing in Eq. (36), for $\omega\simeq 0$ are almost solely determined by the input fields, allowing us to approximate

[TABLE]

In Eq. (37), contrary to first-order case, the term $a^{\mathrm{in}}_{\mathrm{E}\,\omega}$ can give rise to a non-vanishing interference contribution. Focusing only on the terms proportional to $a^{\mathrm{in}}_{\mathrm{I}\,\omega}$, and assuming that $\delta\ll\kappa$,Eq. (37) can be further approximated to give

[TABLE]

From Eqs. (28), (LABEL:eq:4x),(37), it is thus possible to evaluate the output field quadratures (up to second order in the perturbative expansion discussed above). Setting $\lambda=1$ we have

[TABLE]

with

[TABLE]

and analogously for the higher-order mechanical quadrature operators $X_{\omega 1}^{\mathrm{\theta}}$ and $X_{\omega 2}^{\mathrm{\theta}}$.

Considering the thermal input discussed in the main text, and defining the spectrum for the output field as

[TABLE]

with analogous definitions for each perturbative order, the relations given by Eq. (LABEL:eq:9x) allow us to write, up to second order in the detection tone amplitude and for $\omega\simeq 0$, the spectrum of the output noise as

[TABLE]

With the definitions given by Eqs.(28,LABEL:eq:4x, 37), Eq. (42) can be written as

[TABLE]

For $\omega\simeq 0$, we have that the [math]-th order term corresponds to the pure cavity response

[TABLE]

while the terms appearing the third and fourth line of Eq. (43) can be written as

[TABLE]

where we have defined

[TABLE]

The contribution to the output field noise spectrum given by $\kappa_{\mathrm{e}}\braket{X_{-\omega}^{(1)\,\mathrm{\theta}}X_{\omega}^{(1)\,\mathrm{\theta}}}$ can be shown to encode the relevant information about the mechanical quadratures. Again for $\omega\simeq 0$, from (32) we have that

[TABLE]

where

[TABLE]

Appendix C Duan quantity

From the QLEs equations for the mechanical Bogoliubov modes, we evaluate here $\braket{\Delta\bar{X}_{\Sigma}^{2}}+\braket{\Delta\bar{Y}_{\Delta}^{2}}$, which, as we will discuss in the next section, can be shown to correspond to the Duan quantity. From Eqs. (25) (Eqs. (3) of the main text), in the appropriate frame for each mode, we can write the I/O relations for $\beta_{1}$ and $\beta_{2}$ as

[TABLE]

Through Eqs. (4) of the main text, Eqs. (50) can be expressed in terms of mechanical quadrature operators

[TABLE]

and input operators $b_{1\,\mathrm{in}}^{(\dagger)}$, $b_{2\,\mathrm{in}}^{(\dagger)}$, $a_{\mathrm{in}}^{(\dagger)}$ as

[TABLE]

with

[TABLE]

where $\eta^{\theta}_{\pm}=ue^{-i\theta}\mp ve^{i\theta}$. From Eqs. (52) and (53), and assuming $\braket{b^{\dagger}_{1}b_{1}}=\braket{b^{\dagger}_{2}b_{2}}=n_{\mathrm{m}}$, we can evaluate the quadrature variances for the symmetric and antisymmetric modes

[TABLE]

as

[TABLE]

where we have separated the contributions that can be ascribed to the mechanical $\braket{\cdot}_{\mathrm{m}}$ and the optical $\braket{\cdot}_{\mathrm{c}}$ thermal bath. The response integrals can be evaluated analytically, giving (in the limit $\gamma\ll\delta$)

[TABLE]

and

[TABLE]

which, for $\delta\ll\kappa$, allows us to recover the result given in Eq. (14) of the main text.

Appendix D Frequency-shifted quadratures and Duan bound

In order to confirm the presence entanglement between the mechanical resonators, we have to show that the variances of symmetric and antisymmetric quadratures satisfy the Duan bound, which can be expressed in the following form

[TABLE]

where

[TABLE]

with

[TABLE]

where the upper (lower) sign corresponds to $X_{\omega}^{\theta,\Sigma}$ ($X_{\omega}^{\theta,\Delta}$) and

[TABLE]

While the quadratures given in eq. (60) cannot be directly related to the mechanical contribution to the output spectrum, given by $\bar{X}^{\theta}_{\omega,1}$ and $\bar{X}^{\theta}_{\omega,2}$ (see Eq. (9) of the main text), it is possible to show that the following identity holds

[TABLE]

In order to prove Eq. (62), we consider, without loss of generality the case $\theta=0$. The argument of the integral in Eq. (62) can be written in terms of the definitions given in eq. (19) as

[TABLE]

Since the integration has to be performed over the whole frequency domain, upon integration, the frequency in each term can be shifted by the appropriate amount (either $\omega\to\omega+\delta$ or $\omega\to\omega-\delta$), reproducing the result that would be obtained directly evaluating the integral appearing on the rhs of Eq. (62). It is thus clear that the evaluation of the output field spectrum given in Eq. (9) of the main text allows us to determine $S_{\omega}^{\Sigma\,\theta}$, $S_{\omega}^{\Sigma\theta}$, and consequently, the Duan quantity, therefore representing a measure of the degree of entanglement between the mechanical resonators.