Generation and detection of large and robust entanglement between two different mechanical resonators in cavity optomechanics

Entanglement is the distinguishing feature of quantum mechanics and is the physical phenomenon according to which only the properties of the entire system have precise values, while the physical properties of a subsystem can be assigned only in reference to those of the other ones. It is now intensively studied because it corresponds to peculiar nonlocal correlations which allows performing communication and computation tasks with an efficiency which is not achievable classically [1].

Furthermore, for a deeper understanding of the boundary between the classical and quantum world, it is important to investigate up to which macroscopic scale one can observe quantum behavior, and in particular under which conditions entanglement between macroscopic objects, each containing a large number of the constituents, can arise. Entanglement between two atomic ensembles has been successfully demonstrated in [2], while entanglement between two Josephson-junction qubits has been detected in [3, 4]. More recently, macroscopic entanglement has been demonstrated in electro-mechanical systems [5]: continuous variable (CV) entanglement, similar to that considered by Einstein–Podolski and Rosen (EPR) [6], has been generated and detected between the position and momentum of a vibrational mode of a 15 μm diameter Al membrane, and the quadratures of a microwave cavity field, following the theory proposal of [7].

Entanglement between two mechanical resonators (MRs) has been instead demonstrated only at the microscopic level, in the case of two trapped ions [8], and between two single-phonon excitations in nano-diamonds [9]. The realization of this kind of entanglement at the more macroscopic level of micromechanical resonators would be extremely important both for practical and fundamental reasons. In fact, on the one hand, entangled MRs at distant sites could represent an important building block for the implementation of quantum networks for long-distance routing of quantum information [10]; on the other hand, these nonclassical states represent an ideal playground for investigating and comparing decoherence theories and modifications of quantum mechanics at the macroscopic level [11–13].

Many different schemes have been proposed in the literature for entangling two MRs, especially exploiting optomechanical and electromechanical devices [14, 15], in which the two MRs simultaneously interact with one or more electromagnetic cavity fields. References [16–18] considered the steady state of different systems of driven cavities: [16] focused on two mirrors of a ring cavity, while [17] assumed to drive two independent linear cavities with two-mode squeezed light transferring its entanglement to the cavity end-mirrors [18] instead considered a double-cavity scheme in which one cavity couples to the relative motion of two MRs, and the second cavity to their center-of-mass; when the system is appropriately driven by squeezed light, such squeezing is transferred to the two MRs which are then prepared in a stationary EPR-like state. Actually, steady-state entanglement can be achieved, even if at a smaller value, also without squeezed driving, either between two movable mirrors in a Fabry–Perot cavity [19], between two mechanical modes of a single movable mirror [20], or in the case of two semi-transparent membranes interacting with two driven cavity modes [21].

A different approach for generating entangled MRs exploits conditional measurements on light modes entangled or correlated with mechanical degrees of freedom [22–27]. In this case, entanglement is generated at the measurement and it has a finite lifetime which may be severely limited by the interaction of the MRs with their reservoirs. A similar strategy has been provided to enhance the entanglement of two MRs [28]. More recent proposal applied reservoir engineering ideas [29–33] to optomechanical scenarios, by exploiting suitable multi-frequency drivings and optical architectures in order to achieve more robust generation of steady state entanglement between two MRs [34–40], eventually profiting from mechanical nonlinearities and/or parametric driving [41, 42].

In the present paper we propose a novel optomechanical/electromechanical scheme for the generation of remarkably large CV entanglement between two MRs with different frequencies, which is also extremely robust with respect to thermal noise. The scheme is particularly simple, involving only a single, bichromatically-driven, optical cavity mode, and optimally works in a rotating wave approximation (RWA) regime where counter-rotating, non-resonant, terms associated with the bichromatic driving are negligible. The scheme shares some analogies with the reservoir-engineering schemes of [34, 36, 38, 40], but it may be used to generate robust entanglement also in a pulsed regime, in the special case of equal effective couplings at the two sidebands, where the system becomes analogous to the Sørensen–Mølmer scheme for entangling trapped ions in a thermal environment [43]. This latter scheme has been already considered in an optomechanical scenario by Kuzyk et al[44] for entangling dynamically two optical modes via their common interaction with a single MR.

The paper is organized as follows. In section 2 we derive the effective quantum Langevin equations (QLEs) describing the dynamics of the system in the RWA. In section 3 we solve the dynamics in terms of the mechanical Bogoliubov modes of the system [34, 36, 45], derive the steady state of the system in the stable case, and provide simple analytical expressions for the achievable mechanical entanglement, showing its remarkable robustness with respect to temperature. In section 4 we instead consider the special case of equal couplings, when the system can be mapped to the Sørensen–Mølmer scheme [43], in which mechanical entanglement is generated only dynamically and slowly decays to zero at long times. In section 5 we solve and discuss the exact dynamics of the system in order to establish the conditions under which the RWA does not seriously affect the robust generation of large mechanical entanglement. In section 6 we discuss the experimental detection of such entanglement and present some concluding remarks. In the appendices we provide some detail on the dynamical evolution of the system, and present a careful derivation of the linearized QLE in the RWA regime.

As shown in figure 1, we consider an optical cavity mode with resonance frequency ω_c and annihilation operator $\hat{a}$ interacting via the usual optomechanical interaction with two different MRs, with frequencies ω₁ and ω₂ and annihilation operators ${\hat{b}}_{1}$ and ${\hat{b}}_{2}$ respectively. The cavity mode is bichromatically driven at the two frequencies ω₀ + ω₁ and ω₀ − ω₂, with the reference frequency ω₀ detuned from the cavity resonance by a quantity Δ₀ = ω_c − ω₀. If we describe the cavity field in a reference frame rotating at the frequency ω₀, then the system Hamiltonian is given by

$$\begin{eqnarray}\hat{H} & = & {\hslash }{\omega }_{1}{\hat{b}}_{1}^{\dagger }{\hat{b}}_{1}+{\hslash }{\omega }_{2}{\hat{b}}_{2}^{\dagger }{\hat{b}}_{2}+{\hslash }{{\rm{\Delta }}}_{0}{\hat{a}}^{\dagger }\hat{a}\\ & & +{\hslash }\left[{g}_{1}\left({\hat{b}}_{1}+{\hat{b}}_{1}^{\dagger }\right)+{g}_{2}\left({\hat{b}}_{2}+{\hat{b}}_{2}^{\dagger }\right)\right]{\hat{a}}^{\dagger }\hat{a}\\ & & +{\hslash }[({E}_{1}{{\rm{e}}}^{-{\rm{i}}{\omega }_{1}t}+{E}_{2}{{\rm{e}}}^{{\rm{i}}{\omega }_{2}t}){\hat{a}}^{\dagger }+{\rm{h.c.}}].\end{eqnarray} \tag{ 1 }$$

This means that the cavity mode is simultaneously driven on the blue sideband associated with the MR with annihilation operator ${\hat{b}}_{1},$ and on the red sideband associated with the MR with ${\hat{b}}_{2}.$ The nonzero detuning Δ₀ makes the present scheme different from the one studied in the supplementary material of [34] which restricts to the resonant case Δ₀ = 0. Our model is instead related to the scheme proposed by Kuzyk et al [44] for entangling dynamically two optical modes via their common interaction with a single MR: here we will dynamically entangle two MRs via their common interaction with an optical mode.

Zoom In Zoom Out Reset image size

The system dynamics can be efficiently studied by linearizing the optomechanical interaction in the limit of large driving field. In this case the average fields for both cavity, α(t), and mechanical degrees of freedom, β_j(t), are large, and one can simplify the interaction Hamiltonian at lowest order in the field fluctuations

$$\begin{eqnarray}\delta \hat{a}(t) & = & \hat{a}(t)-\alpha (t),\\ \delta {\hat{b}}_{j}(t) & = & {\hat{b}}_{j}(t)-{\beta }_{j}(t).\end{eqnarray} \tag{ 2 }$$

Differently from the typical optomechanical settings in which the steady state average fields are time-independent, here the bichromatic driving induces a time-dependent, periodic steady state average field which, in turn, implies time-dependent effective coupling strengths for the linearized dynamics of the fluctuations. As originally discussed in [46], and detailed in appendix B, approximated dynamical equations for the fluctuation operators $\delta \hat{a}(t)$ and $\delta {\hat{b}}_{j}(t)$ can be derived, in the interaction picture with respect to the Hamiltonian ${\hat{H}}_{0}={\hslash }({\omega }_{1}{\hat{b}}_{1}^{\dagger }{\hat{b}}_{1}+{\omega }_{2}{\hat{b}}_{2}^{\dagger }{\hat{b}}_{2}),$ by neglecting the non-resonant/time-dependent components of the effective linearized interactions. It is possible to prove that this approach is justified when (see equation (B.21))

$$\begin{eqnarray}&&| {g}_{j}\;\displaystyle \frac{{E}_{j}}{{\omega }_{j}}| ,\kappa \ll {\omega }_{j},| {\omega }_{1}-{\omega }_{2}| .\end{eqnarray} \tag{ 3 }$$

The corresponding QLE including thermal noise and dissipation at rates κ and γ_j for the cavity and the mechanical mode $j\in 1,2$ respectively, are

$$\begin{eqnarray}&&\delta \dot{\hat{a}}=-\left(\kappa +{\rm{i}}{\rm{\Delta }}\right)\delta \hat{a}-{\rm{i}}{G}_{1}\delta {\hat{b}}_{1}^{\dagger }-{\rm{i}}{G}_{2}\delta {\hat{b}}_{2}+\sqrt{2\kappa }{\hat{a}}^{{\rm{in}}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\delta {\dot{\hat{b}}}_{1}=-\displaystyle \frac{{\gamma }_{1}}{2}\delta {\hat{b}}_{1}-{\rm{i}}{G}_{1}\delta {\hat{a}}^{\dagger }+\sqrt{{\gamma }_{1}}{\hat{b}}_{1}^{{\rm{in}}},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\delta {\dot{\hat{b}}}_{2}=-\displaystyle \frac{{\gamma }_{2}}{2}\delta {\hat{b}}_{2}-{\rm{i}}{G}_{2}^{*}\delta \hat{a}+\sqrt{{\gamma }_{2}}{\hat{b}}_{2}^{{\rm{in}}},\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}{G}_{1} & = & \displaystyle \frac{{g}_{1}\;{E}_{1}}{{\omega }_{1}-{\rm{\Delta }}+{\rm{i}}\kappa },\\ {G}_{2} & = & -\displaystyle \frac{{g}_{2}\;{E}_{2}}{{\omega }_{2}+{\rm{\Delta }}-{\rm{i}}\kappa },\end{eqnarray} \tag{ 7 }$$

are the (generally complex) linear optomechanical couplings, and ${\hat{a}}^{{\rm{in}}}$ and ${\hat{b}}_{j}^{{\rm{in}}}$ are standard input noise operators with zero mean, whose only nonzero correlation functions are $\left\langle{\hat{a}}^{{\rm{in}}}(t){\hat{a}}^{{\rm{in}}}{({t}^{\prime })}^{\dagger }\right\rangle=\delta (t-{t}^{\prime })$, $\left\langle{\hat{b}}_{j}^{{\rm{in}}}(t){\hat{b}}_{j}^{{\rm{in}}}{({t}^{\prime })}^{\dagger }\right\rangle$ $\;=\;({\bar{n}}_{j}+1)\delta (t-{t}^{\prime })$ and $\left\langle{\hat{b}}_{j}^{{\rm{in}}}{(t)}^{\dagger }{\hat{b}}_{j}^{{\rm{in}}}({t}^{\prime })\right\rangle={\bar{n}}_{j}\delta (t-{t}^{\prime }),$ where ${\bar{n}}_{j}={\left[\mathrm{exp}\left({\hslash }{\omega }_{j}/{k}_{B}T\right)-1\right]}^{-1}$ is the mean thermal phonon number of the jth MR, which we assume to stay at the same environmental temperature T. Moreover, we note that here the new cavity detuning Δ includes the time-independent frequency shift induced by the optomechanical interaction, proportional to the dc component of the average mechanical oscillation amplitude β_j(t), that we here denote with ${\beta }_{j}^{\mathrm{dc}}$ (see appendix B). Specifically

$$\begin{eqnarray}&&{\rm{\Delta }}={{\rm{\Delta }}}_{0}+2\displaystyle \sum _{j=1,2}{g}_{j}{\rm{Re}}[{\beta }_{j}^{{\rm{dc}}}].\end{eqnarray} \tag{ 8 }$$

We will see that the dynamics described by these equations allows to generate large and robust entanglement between the two MRs, either in the steady state or, in a particular parameter regime, during the time evolution with a flat-top pulse driving. We first notice that the system is stable when all the eigenvalues associated with the linearized dynamics of equations (4)–(6) have negative real parts. The stability condition is quite involved in the general case, but it assumes a particularly simple form in the case of equal mechanical dampings, γ₁ = γ₂ = γ. In such a case, the system is stable if and only if

$$\begin{eqnarray}&&{| {G}_{2}| }^{2}\gt {| {G}_{1}| }^{2}-\displaystyle \frac{\kappa \gamma }{2}\left[1+\displaystyle \frac{4{{\rm{\Delta }}}^{2}}{{\left(\gamma +2\kappa \right)}^{2}}\right].\end{eqnarray} \tag{ 9 }$$

This stability condition reduces to the one derived in the supplementary material of [34] in the case Δ = 0. We see that a nonzero detuning generally helps in keeping the system stable.

The coherent dynamics corresponding to the equations (4)–(6), is described by the effective linearized Hamiltonian

$$\begin{eqnarray}{\hat{H}}_{{\rm{eff}}} & = & {\hslash }{\rm{\Delta }}\delta {\hat{a}}^{\dagger }\delta \hat{a}+{\hslash }\left({G}_{1}\delta {\hat{b}}_{1}^{\dagger }+{G}_{2}\delta {\hat{b}}_{2}\right)\delta {\hat{a}}^{\dagger }\\ & & +{\hslash }\left({G}_{1}^{*}\delta {\hat{b}}_{1}+{G}_{2}^{*}\delta {\hat{b}}_{2}^{\dagger }\right)\delta \hat{a}.\end{eqnarray} \tag{ 10 }$$

We can always adjust the phase reference of each MR (which will be determined by a local oscillator which must be used to measure the mechanical quadratures for verifying entanglement) so that we can take both G₁ and G₂ real.

Equation (10) naturally suggests to introduce two effective mechanical modes allowing to simplify the system dynamics. We assume for the moment G₂ > G₁, which is a sufficient condition for stability (see equation (9)), and define

$$\begin{eqnarray}&&{\hat{\beta }}_{1}=\displaystyle \frac{{G}_{2}\delta {\hat{b}}_{1}+{G}_{1}\delta {\hat{b}}_{2}^{\dagger }}{{\mathcal{G}}}=\delta {\hat{b}}_{1}\mathrm{cosh}r+\delta {\hat{b}}_{2}^{\dagger }\mathrm{sinh}r,\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{\hat{\beta }}_{2}=\displaystyle \frac{{G}_{2}\delta {\hat{b}}_{2}+{G}_{1}\delta {\hat{b}}_{1}^{\dagger }}{{\mathcal{G}}}=\delta {\hat{b}}_{2}\mathrm{cosh}r+\delta {\hat{b}}_{1}^{\dagger }\mathrm{sinh}r,\end{eqnarray} \tag{ 12 }$$

where

$$\begin{eqnarray}&&{\mathcal{G}}=\sqrt{{G}_{2}^{2}-{G}_{1}^{2}},\;\;\;\;\mathrm{tanh}r=\displaystyle \frac{{G}_{1}}{{G}_{2}}.\end{eqnarray} \tag{ 13 }$$

Equations (11), (12) define a Bogoliubov unitary transformation of the mechanical mode operators, which can also be written as

$$\begin{eqnarray}{\hat{\beta }}_{\mathrm{1,2}} & = & {{\rm{e}}}^{-r\left(\delta {\hat{b}}_{1}^{\dagger }\delta {\hat{b}}_{2}^{\dagger }-\delta {\hat{b}}_{1}\delta {\hat{b}}_{2}\right)}\delta {\hat{b}}_{\mathrm{1,2}}{{\rm{e}}}^{r\left(\delta {\hat{b}}_{1}^{\dagger }\delta {\hat{b}}_{2}^{\dagger }-\delta {\hat{b}}_{1}\delta {\hat{b}}_{2}\right)}\\ & = & \hat{S}(r)\delta {\hat{b}}_{\mathrm{1,2}}\hat{S}(-r),\end{eqnarray} \tag{ 14 }$$

with $\hat{S}(r)$ the two-mode squeezing operator. The Bogoliubov mode ${\hat{\beta }}_{1}$ describes the 'mechanical dark mode', which does not appear in H_eff, i.e., is decoupled from the cavity mode and therefore is a constant of motion in the absence of damping, while ${\hat{\beta }}_{2}$ is the 'bright' mode interacting with the cavity mode. This is equivalent to say that the dark mode ${\hat{\beta }}_{1}$ is the normal mode of the Hamiltonian dynamics with eigenvalue equal to zero. The other two normal modes of the system will be linear combinations of ${\hat{\beta }}_{2}$ and $\delta \hat{a}.$ The Bogoliubov mode description has been already employed in cavity optomechanics, associated to two optical modes in [34, 45], and to two mechanical modes in [36, 38] (see appendix A for a derivation of the normal modes of the system and a study of its Hamiltonian dynamics).

3.1. Stationary entanglement for different couplings

For a realistic description of the system dynamics we must include cavity decay and mechanical dissipation. It is convenient to rewrite the QLE in terms of the Bogoliubov modes, which in the case when γ₁ = γ₂ ≡ γ assume the simple form

$$\begin{eqnarray}&&\delta \dot{\hat{a}}=-\left(\kappa +{\rm{i}}{\rm{\Delta }}\right)\delta \hat{a}-{\rm{i}}{\mathcal{G}}{\hat{\beta }}_{2}+\sqrt{2\kappa }{\hat{a}}^{{\rm{in}}},\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\dot{\hat{\beta }}}_{1}=-\displaystyle \frac{\gamma }{2}{\hat{\beta }}_{1}+\sqrt{\gamma }{\hat{\beta }}_{1}^{{\rm{in}}},\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{\dot{\hat{\beta }}}_{2}=-\displaystyle \frac{\gamma }{2}{\hat{\beta }}_{2}-{\rm{i}}{\mathcal{G}}\delta \hat{a}+\sqrt{\gamma }{\hat{\beta }}_{2}^{{\rm{in}}},\end{eqnarray} \tag{ 17 }$$

where ${\hat{\beta }}_{j}^{{\rm{in}}},j=1,2,$ are two correlated thermal noise operators whose only nonzero correlation functions are

$$\begin{eqnarray}\left\langle{\hat{\beta }}_{j}^{{\rm{in}}}(t){\hat{\beta }}_{j}^{{\rm{in}}}{({t}^{\prime })}^{\dagger }\right\rangle & = & \left[{\bar{n}}_{j}^{{\rm{eff}}}(r)+1\right]\delta (t-{t}^{\prime }),\\ \left\langle{\hat{\beta }}_{j}^{{\rm{in}}}{(t)}^{\dagger }\;{\hat{\beta }}_{j}^{{\rm{in}}}({t}^{\prime })\right\rangle & = & {\bar{n}}_{j}^{{\rm{eff}}}(r)\delta (t-{t}^{\prime }),\\ \left\langle{\hat{\beta }}_{1}^{{\rm{in}}}(t){\hat{\beta }}_{2}^{{\rm{in}}}({t}^{\prime })\right\rangle & = & \left\langle{\hat{\beta }}_{1}^{{\rm{in}}}{(t)}^{\dagger }\;{\hat{\beta }}_{2}^{{\rm{in}}}{({t}^{\prime })}^{\dagger }\right\rangle=\bar{m}(r)\delta (t-{t}^{\prime }),\end{eqnarray} \tag{ 18 }$$

with the effective mean thermal phonon numbers

$$\begin{eqnarray}&&{\bar{n}}_{1}^{{\rm{eff}}}(r)={\bar{n}}_{1}{\mathrm{cosh}}^{2}r+({\bar{n}}_{2}+1){\mathrm{sinh}}^{2}r,\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{\bar{n}}_{2}^{{\rm{eff}}}(r)={\bar{n}}_{2}{\mathrm{cosh}}^{2}r+({\bar{n}}_{1}+1){\mathrm{sinh}}^{2}r,\end{eqnarray} \tag{ 20 }$$

and the inter-mode correlation

$$\begin{eqnarray}&&\bar{m}(r)=\mathrm{cosh}r\;\mathrm{sinh}r({\bar{n}}_{1}+{\bar{n}}_{2}+1).\end{eqnarray} \tag{ 21 }$$

If ${\gamma }_{1}\ne {\gamma }_{2}$ a dissipative coupling term between the two Bogoliubov modes appears, which however does not have relevant effects because it is proportional to $| {\gamma }_{1}-{\gamma }_{2}|$ which is typically very small with respect to all other damping rates.

The dynamics associated with equations (15)–(17) is simple: the bright mechanical mode ${\hat{\beta }}_{2}$ is cooled by the cavity, while the correlated reservoir create finite correlations between dark and bright modes. In particular, the matrix of correlation for the vector of operators ${\boldsymbol{\beta }}\;=\;$ $({\hat{\beta }}_{1},{\hat{\beta }}_{2},{\hat{\beta }}_{1}^{\dagger },{\hat{\beta }}_{2}^{\dagger }),$ whose elements are ${\{{{\mathcal{C}}}_{\beta }\}}_{j,k}=\left\langle{\{{\boldsymbol{\beta }}\}}_{j}\ {\{{\boldsymbol{\beta }}\}}_{k}\right\rangle$ is given, at the steady state, by

$$\begin{eqnarray}{{\mathcal{C}}}_{\beta }=\left(\begin{array}{cccc}0 & {\bar{m}}_{\beta } & {\bar{n}}_{1}^{{\rm{eff}}}+1 & 0\\ {\bar{m}}_{\beta } & 0 & 0 & {\bar{n}}_{2}^{{\rm{cool}}}+1\\ {\bar{n}}_{1}^{{\rm{eff}}} & 0 & 0 & {\bar{m}}_{\beta }^{*}\\ 0 & {\bar{n}}_{2}^{{\rm{cool}}} & {\bar{m}}_{\beta }^{*} & 0,\end{array}\right)\end{eqnarray} \tag{ 22 }$$

with the number of excitation of the cooled bright mode and the correlations between the two Bogoliubov modes respectively given by

$$\begin{eqnarray}&&{\bar{n}}_{2}^{{\rm{cool}}}(r)={n}_{2}^{{\rm{eff}}}(r)[1-\displaystyle \frac{(1-\epsilon )\;{C}_{-}}{1+{\delta }^{2}+{C}_{-}}],\end{eqnarray} \tag{ 23 }$$

and

$$\begin{eqnarray}&&{\bar{m}}_{\beta }(r)=\bar{m}(r)\displaystyle \frac{2(1+{\rm{i}}\delta )}{2(1+{\rm{i}}\delta )+(1-\epsilon )\;{C}_{-}},\end{eqnarray} \tag{ 24 }$$

where

$$\begin{eqnarray}&&{C}_{-}=\displaystyle \frac{2{{\mathcal{G}}}^{2}}{\gamma \kappa },\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&\epsilon =\displaystyle \frac{\gamma }{\gamma +2\kappa },\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\delta =\displaystyle \frac{2{\rm{\Delta }}}{\gamma +2\kappa },\end{eqnarray} \tag{ 27 }$$

and C₋can be seen as an effective collective optomechanical cooperativity. The steady state correlation matrix can be expressed in terms of the original modes b₁ and b₂ by inverting the Bogoliubov transformation introduced in equations (11) and (12). The result is

$$\begin{eqnarray}{{\mathcal{C}}}_{b} & = & U\ {{\mathcal{C}}}_{\beta }\ {U}^{T}\\ & = & \left(\begin{array}{cccc}0 & {\bar{m}}_{b} & {\bar{n}}_{b1}+1 & 0\\ {\bar{m}}_{b} & 0 & 0 & {\bar{n}}_{b2}+1\\ {\bar{n}}_{b1} & 0 & 0 & {\bar{m}}_{b}\\ 0 & {\bar{n}}_{b2} & {\bar{m}}_{b} & 0,\end{array}\right)\end{eqnarray} \tag{ 28 }$$

with

$$\begin{eqnarray}U=\left(\begin{array}{cccc}\mathrm{cosh}r & 0 & 0 & -\mathrm{sinh}r\\ 0 & \mathrm{cosh}r & -\mathrm{sinh}r & 0\\ 0 & -\mathrm{sinh}r & \mathrm{cosh}r & 0\\ -\mathrm{sinh}r & 0 & 0 & \mathrm{cosh}r,\end{array}\right)\end{eqnarray} \tag{ 29 }$$

and where now

$$\begin{eqnarray}{\bar{n}}_{b1} & = & {\bar{n}}_{1}^{{\rm{eff}}}+{\mathrm{sinh}}^{2}r(1+{\bar{n}}_{1}^{{\rm{eff}}}+{\bar{n}}_{2}^{{\rm{cool}}})\\ & & -2\mathrm{cosh}r\ \mathrm{sinh}r\ {\rm{Re}}({\bar{m}}_{\beta }),\\ {\bar{n}}_{b2} & = & {\bar{n}}_{2}^{{\rm{cool}}}+{\mathrm{sinh}}^{2}r(1+{\bar{n}}_{1}^{{\rm{eff}}}+{\bar{n}}_{2}^{{\rm{cool}}})\\ & & -2\mathrm{cosh}r\ \mathrm{sinh}r\ {\rm{Re}}({\bar{m}}_{\beta }),\\ {\bar{m}}_{b} & = & {\mathrm{cosh}}^{2}r\ {\bar{m}}_{\beta }+{\mathrm{sinh}}^{2}r\ {\bar{m}}_{\beta }^{*}\\ & & -\mathrm{cosh}r\ \mathrm{sinh}r(1+{\bar{n}}_{1}^{{\rm{eff}}}+{\bar{n}}_{2}^{{\rm{cool}}}).\end{eqnarray} \tag{ 30 }$$

The entanglement between modes b₁ and b₂, measured by means of the logarithmic negativity [47, 48], can be easily expressed in terms of these matrix elements as [49]

$$\begin{eqnarray}{E}_{N} & = & {\rm{max}}[0,-\mathrm{ln}\nu ],\\ \nu & = & 1+{\bar{n}}_{b1}+{\bar{n}}_{b2}-\sqrt{4{| {\bar{m}}_{b}| }^{2}+{({\bar{n}}_{b1}-{\bar{n}}_{b2})}^{2}}.\end{eqnarray} \tag{ 31 }$$

When the collective cooperativity C₋is sufficiently large, i.e., ${C}_{-}\gg \bar{m}(r),$ then ${\bar{m}}_{\beta }(r)$ is negligible (see equation (24)). This is the working regime in which we are particularly interested, because in this case, the second Bogoliubov mode can be cooled close to its ground state (${\bar{n}}_{2}^{{\rm{cool}}}(r)\ll {n}_{2}^{{\rm{eff}}}(r)$), corresponding to an entangled state for the original mechanical modes. In this case the steady state correlation matrix for the Bogoliubov modes, in equation (22), reduces to the correlation matrix of a state given by the product of two thermal states with occupancies ${\bar{n}}_{1}^{{\rm{eff}}}(r)$ and ${\bar{n}}_{2}^{{\rm{cool}}}(r)$ respectively. For the two MR of interest, associated with the operator ${\hat{b}}_{1}$ and ${\hat{b}}_{2},$ such a state is just a two-mode squeezed thermal state [50]

$$\begin{eqnarray}&&{\hat{\rho }}_{\mathrm{1,2}}=\hat{S}(r){\hat{\rho }}_{{\bar{n}}_{1}^{{\rm{eff}}}(r),{\rm{th}}}\otimes {\hat{\rho }}_{{\bar{n}}_{2}^{{\rm{cool}}}(r),{\rm{th}}}\hat{S}(-r),\end{eqnarray} \tag{ 32 }$$

where $\hat{S}(r)$ is given in equation (14), and

$$\begin{eqnarray}&&{\hat{\rho }}_{\bar{n},{\rm{th}}}=\displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{{\bar{n}}^{n}}{{(1+\bar{n})}^{n+1}}| n\rangle \langle n| \end{eqnarray} \tag{ 33 }$$

is the density matrix of the thermal equilibrium state of a resonator with occupancy $\bar{n}.$ Such a state is entangled for sufficiently large r and not too large mean thermal excitation number.

This prediction of large stationary entanglement is confirmed in figure 2, where we plot the time evolution of the entanglement between the two MRs, quantified in terms of the logarithmic negativity E_N, obtained from the solution of equations (4)–(6). Figure 2 refers to an experimentally achievable set of parameters, γ = 10 s⁻¹, κ = 10⁵ s⁻¹, G₂ = 10⁵ s⁻¹, Δ = 10³ s⁻¹, and to different values of mean thermal phonon numbers ${\bar{n}}_{1},{\bar{n}}_{2},$ and of the ratio G₁/G₂. We see that remarkable values of E_N are achieved at low temperatures, and that stationary mechanical entanglement is quite robust with respect to temperature because one has an appreciable value of E_N ≃ 0.32 even for ${\bar{n}}_{1}=2000,{\bar{n}}_{2}=1000$. The time to reach the steady state is essentially given by the inverse of the cooling rate of the bright Bogoliubov mode, which is approximately given by t_s ≃ (κ² + Δ²)/(${{\mathcal{G}}}^{2}\kappa$) (see equations (15)–(17)).

Zoom In Zoom Out Reset image size

Equation (32) suggests that one could achieve large stationary entanglement between the two MRs by taking a large two-mode squeezing parameter r, and a large collective cooperativity C₋$\gg$ 1 in order to significantly cool the bright Bogoliubov mode. However the corresponding optimization of the system parameters, and especially of the two couplings G₁ and G₂, is far from being trivial. In fact, r increases when ${G}_{1}\to {G}_{2},$ which however implies, at a fixed value of G₂, a decreasing value of ${\mathcal{G}}$ and therefore of C₋(see equations (13) and (25)); moreover increasing r has also the unwanted effect of increasing ${\bar{m}}_{\beta }(r)$ that is the correlations between the two Bogoliubov modes (see equations (21) and (24)).

However, a judicious choice of parameters is possible, allowing to get very large stationary mechanical entanglement, even in the presence of non-negligible values of the thermal occupancies ${\bar{n}}_{1}$ and ${\bar{n}}_{2}.$ At a given value of G₁, this is obtained by taking a sufficiently large value of the associated single-mode cooperativity, ${C}_{1}=2{G}_{1}^{2}/\kappa \gamma \gg 1,$ and correspondingly optimizing the value of G₂, i.e., of r. In fact, the logarithmic negativity associated with the stationary state of equation (32) can be evaluated in terms of the parameter

$$\begin{eqnarray}{\nu | }_{{\bar{m}}_{\beta }\to 0} & = & [{\bar{n}}_{+}(r)+1]({\mathrm{cosh}}^{2}r+{\mathrm{sinh}}^{2}r)\\ & & -\sqrt{{\bar{n}}_{-}{(r)}^{2}+4{[{\bar{n}}_{+}(r)+1]}^{2}{\mathrm{sinh}}^{2}r{\mathrm{cosh}}^{2}r},\end{eqnarray} \tag{ 34 }$$

where ${\bar{n}}_{\pm }(r)={\bar{n}}_{1}^{{\rm{eff}}}(r)\pm {\bar{n}}_{2}^{{\rm{cool}}}(r),$ and ${\bar{n}}_{2}^{{\rm{cool}}}(r)$ can be explicitly rewritten in terms of the cooperativity C₁ as

$$\begin{eqnarray}{\bar{n}}_{2}^{{\rm{cool}}}(r) & = & \left[{\bar{n}}_{2}{\mathrm{cosh}}^{2}r+({\bar{n}}_{1}+1){\mathrm{sinh}}^{2}r\right]\\ & & \times [1-\displaystyle \frac{(1-\epsilon )\;{C}_{1}}{{\mathrm{sinh}}^{2}r(1+{\delta }^{2})+{C}_{1}}].\end{eqnarray} \tag{ 35 }$$

The dependence of E_N versus r, for given values of ${C}_{1}{\bar{n}}_{1}$ and ${\bar{n}}_{2},$ shows a maximum and then decays to zero for large r (see figure 3 which refers to C₁ = 2 × 10⁴ and ${\bar{n}}_{1}=200,{\bar{n}}_{2}=100$). This behavior is described by a very simple approximated expression valid in the limit ${C}_{1}\gg {{\rm{e}}}^{2r}\gg {{\rm{e}}}^{-2r},$ with not very large ${\bar{n}}_{\mathrm{1,2}},$ and when δ, → 0 (corresponding to $\gamma ,{\rm{\Delta }}\ll \kappa$ ),

$$\begin{eqnarray}&&\nu \sim 2{{\rm{e}}}^{-2r}+(1+{\bar{n}}_{1}+{\bar{n}}_{2})\displaystyle \frac{{{\rm{e}}}^{2r}}{4{C}_{1}}\end{eqnarray} \tag{ 36 }$$

which exhibits a minimum (hence corresponding to maximum entanglement) as a function of r at

$$\begin{eqnarray}&&r\simeq {r}^{{\rm{opt}}}=\displaystyle \frac{1}{4}\mathrm{ln}\left(\displaystyle \frac{8{C}_{1}}{{\bar{n}}_{1}+{\bar{n}}_{2}+1}\right),\end{eqnarray} \tag{ 37 }$$

given by ${\nu }^{{\rm{opt}}}\sim \sqrt{\frac{2(1+{\bar{n}}_{1}+{\bar{n}}_{2})}{{C}_{1}}}.$ For values of r much larger or much smaller than this value, the resonators may not be entangled. When r is increased to very large values $r\gg {r}^{{\rm{opt}}},{\mathcal{G}}$ is reduced and the cooling dynamics becomes slow as compared to the standard mechanical dissipation, which takes place at rate $\sim \gamma \ ({\bar{n}}_{j}+1),$ so that the correlations between the MRs cannot be efficiently generated. On the other hand, at small $r\ll {r}^{{\rm{opt}}}$ the Bogoliubov modes are essentially equal to the original modes, so that the cavity cools only the second resonator, and also in this case mechanical entanglement can not be observed. Figure 3 also shows that the simplified expression of equation (36) provides a simple but valid approximation for large C₁ and a very good estimate of the optimal value of the two-mode squeezing parameter r, i.e., of G₁/G₂, given by equation (37). The corresponding value of the logarithmic negativity is

$$\begin{eqnarray}&&{E}_{N}\sim \displaystyle \frac{1}{2}\mathrm{ln}\left[\displaystyle \frac{{C}_{1}}{2(1+{\bar{n}}_{1}+{\bar{n}}_{2})}\right]\end{eqnarray} \tag{ 38 }$$

and shows that once that the ratio G₁/G₂ is optimized, the achievable stationary entanglement between the two MRs increases with increasing ${C}_{1}/({\bar{n}}_{1}+{\bar{n}}_{2}).$

Zoom In Zoom Out Reset image size

The above analysis of the stationary entanglement of the two MRs extends the results of [34] in various directions. First of all, our model extends to the case of nonzero detuning Δ a model discussed in the Supplementary material of [34]. We see that a nonzero detuning has a limited effect of the dynamic of entanglement generation, providing only an effective increase of ${\bar{n}}_{2}^{{\rm{cool}}},$ which however becomes negligible as soon as ${\rm{\Delta }}\ll \kappa$ (see equation (35)). Moreover, [34] provided an explicit expression for E_N only for the case of negligible thermal occupancies and not too large values of r, while the present discussion applies for arbitrary values of $r,{\bar{n}}_{1}$ and ${\bar{n}}_{2}.$

In the special case of equal couplings G₁ = G₂ ≡ G, i.e., ${\mathcal{G}}=0,$ the Bogoliubov modes cannot be defined anymore and the description of the preceding section cannot be applied. The dynamics is nonetheless interesting and still allows for the generation of appreciable entanglement between the two MRs, even though only at finite times and not in the stationary state. We notice that in this special case, our scheme becomes analogous to that of [44], that showed that two appropriately driven optical modes can be entangled with a pulsed scheme by their common interaction with a MR. More precisely, the QLE of equations (4)–(6) are the same as those studied in [44] but now referred to two mechanical modes coupled to the same optical mode, i.e., with exchanged roles between optical and mechanical degrees of freedom.

The physical mechanism at the basis of the generation of dynamical entanglement can be understood by looking at the Hamiltonian evolution of the system at equal couplings. Such mechanism essentially coincides with the one proposed for entangling internal states of trapped ions by Milburn [51] and by Sørensen and Mølmer [43], and first applied to an optomechanical setup by Kuzyk et al [44]. In the present case, the common interaction with the bichromatically driven optical mode dynamically entangles the two MRs, and at special values of the interaction time the optical mode is decoupled from the two MRs and mechanical entanglement can be strong.

At equal couplings it is convenient to rewrite the effective Hamiltonian after linearization of equation (10) in terms of mechanical and optical quadratures, using the expressions $\delta {\hat{b}}_{j}=({\hat{x}}_{j}+{\rm{i}}{\hat{p}}_{j})/\sqrt{2}j=1,2,$ and $\delta \hat{a}=(\hat{X}+{\rm{i}}\hat{Y})/\sqrt{2}.$ One gets

$$\begin{eqnarray}&&{\hat{H}}_{{\rm{eff}}}=\displaystyle \frac{{\hslash }{\rm{\Delta }}}{2}\left({\hat{X}}^{2}+{\hat{Y}}^{2}\right)+{\hslash }G\sqrt{2}\left({\hat{x}}_{+}\hat{X}-{\hat{p}}_{-}\hat{Y}\right),\end{eqnarray} \tag{ 39 }$$

where ${\hat{x}}_{\pm }=({\hat{x}}_{1}\pm {\hat{x}}_{2})/\sqrt{2},{\hat{p}}_{\pm }=({\hat{p}}_{1}\pm {\hat{p}}_{2})/\sqrt{2}$ are linear combinations of the two position and momentum operators of the two MRs. The Heisenberg evolution of these latter mechanical operators can be solved in a straightforward way, by exploiting the fact that ${\hat{x}}_{+}$ and ${\hat{p}}_{-}$ are two commuting conserved observables. One gets (see also [43, 44, 51])

$$\begin{eqnarray}{\hat{x}}_{+}(t) & = & {\hat{x}}_{+}(0),\;\;\;\;{\hat{p}}_{-}(t)={\hat{p}}_{-}(0),\\ {\hat{x}}_{-}(t) & = & {\hat{x}}_{-}(0)+\displaystyle \frac{2{G}^{2}}{{{\rm{\Delta }}}^{2}}\left(\mathrm{sin}{\rm{\Delta }}t-{\rm{\Delta }}t\right){\hat{p}}_{-}(0)\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&+\displaystyle \frac{2{G}^{2}}{{{\rm{\Delta }}}^{2}}\left(1-\mathrm{cos}{\rm{\Delta }}t\right){\hat{x}}_{+}(0)\\ &&-\displaystyle \frac{G\sqrt{2}}{{\rm{\Delta }}}\mathrm{sin}{\rm{\Delta }}t\hat{Y}(0)+\displaystyle \frac{G\sqrt{2}}{{\rm{\Delta }}}\left(1-\mathrm{cos}{\rm{\Delta }}t\right)\hat{X}(0),\\ &&{\hat{p}}_{+}(t)={\hat{p}}_{+}(0)-\displaystyle \frac{2{G}^{2}}{{{\rm{\Delta }}}^{2}}\left(\mathrm{sin}{\rm{\Delta }}t-{\rm{\Delta }}t\right){\hat{x}}_{+}(0)\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}&&+\displaystyle \frac{2{G}^{2}}{{{\rm{\Delta }}}^{2}}\left(1-\mathrm{cos}{\rm{\Delta }}t\right){\hat{p}}_{-}(0)\\ &&-\displaystyle \frac{G\sqrt{2}}{{\rm{\Delta }}}\mathrm{sin}{\rm{\Delta }}t\hat{X}(0)-\displaystyle \frac{G\sqrt{2}}{{\rm{\Delta }}}\left(1-\mathrm{cos}{\rm{\Delta }}t\right)\hat{Y}(0).\end{eqnarray} \tag{ 42 }$$

Relevant interaction times are those when the MR dynamics decouple from that of the optical cavity, and this occurs at t_m = 2m π/Δm = 1, 2, ..., where

$$\begin{eqnarray}&&{\hat{x}}_{-}({t}_{m})={\hat{x}}_{-}(0)-2m\pi \displaystyle \frac{2{G}^{2}}{{{\rm{\Delta }}}^{2}}{\hat{p}}_{-}(0),\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&{\hat{p}}_{+}(t)={\hat{p}}_{+}(0)+2m\pi \displaystyle \frac{2{G}^{2}}{{{\rm{\Delta }}}^{2}}{\hat{x}}_{+}(0).\end{eqnarray} \tag{ 44 }$$

This map describes a stroboscopic evolution in which the two MRs become more and more entangled, because it corresponds to the application of the unitary operator

$$\begin{eqnarray}{U}_{m} & = & \mathrm{exp}\left[-{\rm{i}}\displaystyle \frac{2\pi {G}^{2}m}{{{\rm{\Delta }}}^{2}}\left({\hat{x}}_{+}^{2}+{\hat{p}}_{-}^{2}\right)\right]\\ & = & \mathrm{exp}\left[-{\rm{i}}\displaystyle \frac{2\pi {G}^{2}m}{{{\rm{\Delta }}}^{2}}\left(\delta {\hat{b}}_{1}^{\dagger }\delta {\hat{b}}_{1}+\delta {\hat{b}}_{2}^{\dagger }\delta {\hat{b}}_{2}+1+\delta {\hat{b}}_{1}^{\dagger }\delta {\hat{b}}_{2}^{\dagger }+\delta {\hat{b}}_{1}\delta {\hat{b}}_{2}\right)\right].\end{eqnarray} \tag{ 45 }$$

This ideal behavior is significantly modified by the inclusion of damping and noise, especially the one associated with the cavity mode, which acts on the faster timescale 1/κ and seriously affects the cavity-mediated interaction between the two MRs, as soon as κ becomes comparable to Δ. Mechanical entanglement is large for large G/Δ and we expect well distinct peaks for E_N at interaction times t_m, in the ideal parameter regime $G\gg {\rm{\Delta }}\gg \kappa .$ In the more realistic regime in which GΔ and κ are comparable, the peaks will be washed out, but we still expect an appreciable value for the mechanical entanglement for a large interval of interaction times. This is confirmed by the numerical solution of the time evolution associated with the QLE shown in figure 4, which refers to the parameter set γ = 10 s⁻¹, κ = 10⁵ s⁻¹, G = 10⁵${{\rm{s}}}^{-1},{\bar{n}}_{1}=200,{\bar{n}}_{2}=100,$ and to three different values of the detuning, Δ = 10³ s⁻¹ (black dashed line), Δ = 10⁴ s⁻¹ (red full line), and Δ = 10⁵ s⁻¹ (blue full line). We see that an appreciable value of E_N (even though smaller than the one achievable at the same ${\bar{n}}_{1}$ and ${\bar{n}}_{2}$ after the optimization of G₁/G₂ of the previous Section) is reached for a large interval of interaction times t. Therefore even at equal couplings (and nonzero detuning) one can entangle the two resonators with a pulsed experiment. Mechanical entanglement instead vanishes in the stationary state.

Zoom In Zoom Out Reset image size

The derivation of the effective linearized dynamics of appendix B suggests that the counter-rotating terms that we have neglected may play an important role when the mechanical frequencies are not too large with respect to the other parameters (see also the comments in the supplementary material of [34]). It is therefore interesting to study their effect by comparing the above predictions, both in the case of G₂ > G₁ and in the case of equal couplings, to the solution of the exact QLE obtained without neglecting the various time-dependent terms.

In appendix B we describe the derivation of the effective linearized equations that we have studied in the preceding sections and that is based on the elimination of fast rotating terms and on the expansion of the linearized coupling strength at lowest order in g_j. Here we analyze the limit of validity of these approximations by solving numerically the system dynamics with the inclusion of the non-resonant terms expanded at different orders in powers of g_j. In figures 5 and 6 the red lines are evaluated without the non-resonant terms (i.e., the treatment of the preceding sections), while the green and the blue ones take into account the full dynamics. In particular the green lines are computed by expanding the average fields α(t) and β_j(t) (that have been introduced in equation (2) and discussed in appendix B), at the lowest relevant order in powers of g, while for the blue ones they have been expanded at sixth order in powers of g. Moreover, the green line results are found considering only the steady state solution for α(t) and β_j(t), while the blue lines are computed taking into account their full dynamics (that includes also the transient regime before the steady state is reached) with initial condition α(0) = β_j(0) = 0.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figure 5 we compare the time evolution of the entanglement evaluated with and without the time-dependent terms when G₂ > G₁. The parameters used in these plots are consistent with those used in figure 2. Specifically the three red curves in figures 5 (a), (b) and (c), that are barely visible because almost entirely covered by the green curves, are equal to the three lowest curves in figure 2. We observe that the green and the red lines are always very close, meaning that the linearized RWA treatment is a very good approximation of the full dynamics when α(t) and β_j(t) can be expanded at lowest order in g. Nevertheless, we note that if the mechanical frequencies are not large enough and higher order terms are taken into account together with the full dynamics of α(t) and β_j(t), then the results can be significantly different as described by the blue curves. Specifically, the solid-blue lines are evaluated for sufficiently large values of the mechanical frequencies so that the condition in equation (3) is well fulfilled, and the effective linearized RWA dynamics recovers with significant accuracy the one determined with the inclusion of the non-resonant terms The dashed-blue lines are instead evaluated for smaller frequencies. In this case it is evident that the non-resonant terms have a significant role in the system dynamics and that the lowest order expansion of the coefficients α(t) and β_j(t) does not provide an accurate description. We note that according to equation (3), in order to eliminate the fast rotating terms, the ratios ω_j/G_j have to be much larger than one. Although the dashed-blue curves are evaluated for a ratio ω₁/G₁ of roughly 25, which can be considered significantly large, we have found, indeed, that it is not enough for a faithful approximation of the system dynamics with the model discussed in the preceding sections. The conclusive analysis of these cases would, possibly, require a non-perturbative approach that is beyond the scope of the present work. A final remark is in order. We have verified that the discrepancy between the dashed-blue lines and the red ones is due to the combined effect of the higher order terms and of the transient initial dynamics of α(t) and β_j(t). Specifically, when we consider either the lowest order terms and the transient dynamics, or the higher order terms and only the steady state of α(t) and β_j(t), the corresponding results for the entanglement dynamics are very similar to the red lines.

In figure 6 we study the case of equal couplings G₁ = G₂. In this case solid and dashed lines differ in the values of the cavity detuning Δ. In general larger Δ (dashed lines) corresponds to smaller entanglement, and the results evaluated by including the counter-rotating terms tends to exhibit larger entanglement than the corresponding ones obtained without the non-resonant terms. The solid curves are found with smaller Δ. In this case red, green and blue lines are very close when the mechanical dissipation is sufficiently large as in figure 6 (a). Larger discrepancies are found when the mechanical dissipation is reduced as in figures 6 (b) and (c), especially at relatively large time. We observe in fact that, while the red curves for the entanglement decay to zero at large time, the corresponding green and blue lines seem to approach a finite sizable value. As shown by the insets, when this different behavior is observed, the average photon number in the cavity tends to diverge. This is a signature of the fact that the full dynamics including counter-rotating terms is actually unstable, even though the RWA dynamics without these terms is stable (see equation (9)). We have confirmed the unstable nature of the time-dependent dynamics by calculating the Floquet exponents of the dynamical equations of the system. In fact, when α(t) and β_j(t) are considered in their steady state, one has a system of linear differential equations with periodic, time-dependent coefficients (see appendix B), and the Floquet theory can be applied in this case [52]; we have verified that for the parameters of figure 6 there is always at least one positive Floquet exponent, meaning that the system is unstable. This implies that, in general, the corresponding results are well-grounded only for relatively short time until the populations are not exceedingly large. On the other hand, our results show that in a pulsed experiment with the parameters of figure 6, these instabilities do not constitute a serious hindrance to the creation of significant entanglement at finite times.

Therefore, when the mechanical frequencies are sufficiently large (${\omega }_{j}\gtrsim {10}^{2}\kappa$) (and, limited only to the case of equal couplings, when also mechanical damping is not too small), the effective linearized RWA dynamics obtained by neglecting the counter-rotating terms approximates with very good accuracy the full system dynamics.

We finally discuss how to detect the generated mechanical entanglement between the two MRs at different frequencies. The present entanglement describes EPR-like correlations between the quadratures of the two MRs and therefore we need to perform homodyne-like detection of these quadratures. In the linearized regime we are considering the state of the two MRs is a Gaussian CV state, which is fully characterized by the matrix of all second-order correlations between the mechanical quadratures. Therefore from the measurement of these correlations one can extract the logarithmic negativity E_N. One does not typically have direct access to the mechanical quadratures, but one can exploit the currently available possibility to perform low-noise and highly efficient homodyne detection of optical and microwave fields, and implement an efficient transfer of the mechanical phase-space quadratures onto the optical/microwave field.

As suggested in [53] and then implemented in the electromechanical entanglement experiment of [5], the motional quadratures of a MR can be read by homodyning the output of an additional 'probe' cavity mode. In particular, if the readout cavity mode is driven by a much weaker laser so that its back-action on the mechanical mode can be neglected, and resonant with the first red sideband of the mode, i.e., with a detuning ${{\rm{\Delta }}}_{j}^{p}={\omega }_{j},j=1,2,$ the probe mode adiabatically follows the MR dynamics, and the output of the readout cavity ${a}_{j}^{\mathrm{out}}$ is given by (see figure 1) [53]

$$\begin{eqnarray}&&{a}_{j}^{\mathrm{out}}={\rm{i}}\displaystyle \frac{{G}_{j}^{p}}{\sqrt{\kappa }}\delta {b}_{j}+{a}_{j}^{\mathrm{in}},\;\;\;j=1,2,\end{eqnarray} \tag{ 46 }$$

with ${G}_{j}^{p}$ the very small optomechanical coupling with the probe mode. Therefore using a probe mode for each MR, changing the phases of the corresponding local oscillator, and measuring the correlations between the probe mode outputs, one can then detect all the entries of the correlation matrix and from them numerically extract the logarithmic negativity E_N.

6.1. Concluding remarks

This work has been supported by the European Commission (ITN-Marie Curie project cQOM, Grant No. 290161, and FET-Open Project iQUOEMS, Grant No. 323924), by MIUR (PRIN 2011).

It is straightforward to see that the diagonal form of the interaction Hamiltonian of equation (10) is

$$\begin{eqnarray}&&{\hat{H}}_{{\rm{eff}}}={\hslash }{\lambda }_{0}{\hat{\beta }}_{1}^{\dagger }{\hat{\beta }}_{1}+{\hslash }{\lambda }_{1}{\hat{\alpha }}_{1}^{\dagger }{\hat{\alpha }}_{1}+{\hslash }{\lambda }_{2}{\hat{\alpha }}_{2}^{\dagger }{\hat{\alpha }}_{2},\end{eqnarray} \tag{ A.1 }$$

where

$$\begin{eqnarray}&&{\hat{\alpha }}_{1}=\mathrm{cos}\theta {\hat{\beta }}_{2}+\mathrm{sin}\theta \delta \hat{a},\end{eqnarray} \tag{ A.2 }$$

$$\begin{eqnarray}&&{\hat{\alpha }}_{2}=\mathrm{cos}\theta \delta \hat{a}-\mathrm{sin}\theta {\hat{\beta }}_{2},\end{eqnarray} \tag{ A.3 }$$

define the other two normal modes together with the dark mode ${\hat{\beta }}_{1},$ introduced in equation (11), with θ defined by the condition $\mathrm{tan}2\theta =-2{\mathcal{G}}/{\rm{\Delta }},$ while the eigenvalues are given by ${\lambda }_{0}=0{\lambda }_{1}\;=\;$ $\left({\rm{\Delta }}-\tilde{{\rm{\Delta }}}\right)/2{\lambda }_{2}\;=\;$ $\left({\rm{\Delta }}+\tilde{{\rm{\Delta }}}\right)/2,$ with $\tilde{{\rm{\Delta }}}=\sqrt{{{\rm{\Delta }}}^{2}+4{{\mathcal{G}}}^{2}}.$

The normal modes allows to understand the dynamics in the absence of optical and mechanical damping processes. In fact, from equation (A.1) one can easily derive the Heisenberg evolution of the mechanical bosonic operators. By inverting equations (11), (12) one has, $\delta {\hat{b}}_{1}(t)=\mathrm{cosh}r{\hat{\beta }}_{1}(t)-\mathrm{sinh}r{\hat{\beta }}_{2}^{\dagger }(t)\delta {\hat{b}}_{2}(t)\;=\;$ $\mathrm{cosh}r{\hat{\beta }}_{2}(t)-\mathrm{sinh}r{\hat{\beta }}_{1}^{\dagger }(t)$ and using ${\alpha }_{j}(t)\;=\;$ ${{\rm{e}}}^{{\rm{i}}\lambda t}{\alpha }_{j}(0),j=1,2,$ and β₁(t) = β₁(0), one gets

$$\begin{eqnarray}\delta {\hat{b}}_{1}(t) & = & \mathrm{cosh}r{\hat{\beta }}_{1}(0)\\ & & -\mathrm{sinh}r\mathrm{exp}\left[{\rm{i}}\displaystyle \frac{{\rm{\Delta }}t}{2}\right]\left[\mathrm{cos}\displaystyle \frac{\tilde{{\rm{\Delta }}}t}{2}-{\rm{i}}\;\mathrm{cos}\;2\theta \mathrm{sin}\displaystyle \frac{\tilde{{\rm{\Delta }}}t}{2}\right]{\hat{\beta }}_{2}^{\dagger }(0)\\ & & -{\rm{i}}\;\mathrm{sinh}r\mathrm{sin}2\theta \mathrm{exp}\left[{\rm{i}}\displaystyle \frac{{\rm{\Delta }}t}{2}\right]\mathrm{sin}\displaystyle \frac{\tilde{{\rm{\Delta }}}t}{2}\delta \hat{a}(0)\end{eqnarray} \tag{ A.4 }$$

$$\begin{eqnarray}\delta {\hat{b}}_{2}(t) & = & -\mathrm{sinh}r{\hat{\beta }}_{1}^{\dagger }(0)\\ & & +\mathrm{cosh}r\mathrm{exp}\left[-{\rm{i}}\displaystyle \frac{{\rm{\Delta }}t}{2}\right]\left[\mathrm{cos}\displaystyle \frac{\tilde{{\rm{\Delta }}}t}{2}+{\rm{i}}\;\mathrm{cos}\;2\theta \mathrm{sin}\displaystyle \frac{\tilde{{\rm{\Delta }}}t}{2}\right]{\hat{\beta }}_{2}(0)\\ & & +{\rm{i}}\;\mathrm{cosh}r\mathrm{sin}2\theta \mathrm{exp}\left[-{\rm{i}}\displaystyle \frac{{\rm{\Delta }}t}{2}\right]\mathrm{sin}\displaystyle \frac{\tilde{{\rm{\Delta }}}t}{2}\delta \hat{a}(0).\end{eqnarray} \tag{ A.5 }$$

We now look for special time instants at which the two mechanical modes can be strongly entangled. A necessary condition for such dynamical entanglement is that at these times, the cavity mode must be decoupled from the mechanical modes and equations (A.4), (A.5) show that it occurs when $\mathrm{sin}\tilde{{\rm{\Delta }}}t/2=0,$ i.e., ${t}_{p}=2\pi p/\sqrt{{{\rm{\Delta }}}^{2}+4{{\mathcal{G}}}^{2}}p=1,2,\ldots .$ At these time instants one has

$$\begin{eqnarray}\delta {\hat{b}}_{1}({t}_{p}) & = & \left[{\mathrm{cosh}}^{2}r-{{\rm{e}}}^{{\rm{i}}{\phi }_{p}}{\mathrm{sinh}}^{2}r\right]\delta {\hat{b}}_{1}(0)\\ & & +\mathrm{sinh}r\mathrm{cosh}r\left(1-{{\rm{e}}}^{{\rm{i}}{\phi }_{p}}\right)\delta {\hat{b}}_{2}^{\dagger }(0),\end{eqnarray} \tag{ A.6 }$$

$$\begin{eqnarray}\delta {\hat{b}}_{2}^{\dagger }({t}_{p}) & = & \left[{{\rm{e}}}^{{\rm{i}}{\phi }_{p}}{\mathrm{cosh}}^{2}r-{\mathrm{sinh}}^{2}r\right]\delta {\hat{b}}_{2}^{\dagger }(0)\\ & & -\mathrm{sinh}r\mathrm{cosh}r\left(1-{{\rm{e}}}^{{\rm{i}}{\phi }_{p}}\right)\delta {\hat{b}}_{1}(0),\end{eqnarray} \tag{ A.7 }$$

where ${\phi }_{p}=\pi p\left(1+{\rm{\Delta }}/\tilde{{\rm{\Delta }}}\right).$ In particular, if ${{\rm{e}}}^{{\rm{i}}{\phi }_{p}}=-1$ one gets

$$\begin{eqnarray}&&\delta {\hat{b}}_{1}({t}_{p})=\mathrm{cosh}2r\delta {\hat{b}}_{1}(0)+\mathrm{sinh}2r\delta {\hat{b}}_{2}^{\dagger }(0),\end{eqnarray} \tag{ A.8 }$$

$$\begin{eqnarray}&&\delta {\hat{b}}_{2}^{\dagger }({t}_{p})=-\mathrm{cosh}2r\delta {\hat{b}}_{2}^{\dagger }(0)-\mathrm{sinh}2r\delta {\hat{b}}_{1}(0),\end{eqnarray} \tag{ A.9 }$$

i.e., the state of the two MRs at time t_p is the result of the application of the two-mode squeezing operator with squeezing parameter $2r,\hat{S}(2r)$ (see equation (14)) to their initial state. In the usual case of an initial thermal state for the two MRs with mean thermal phonon numbers ${\bar{n}}_{j},$ the state at time t_p is therefore a two-mode squeezed thermal state [50] (see equation (32)), with logarithmic negativity [47, 48]

$$\begin{eqnarray}{E}_{N}({t}_{p}) & = & -\displaystyle \frac{1}{2}\mathrm{ln}\left[{0ex}{3ex},{\bar{n}}_{-}^{2}+{({\bar{n}}_{+}+1)}^{2}\mathrm{cosh}8r\right.\\ & & -\left.\sqrt{{({\bar{n}}_{+}+1)}^{4}{\mathrm{sinh}}^{2}8r+4{\bar{n}}_{-}^{2}{({\bar{n}}_{+}+1)}^{2}{\mathrm{cosh}}^{2}4r}{0ex}{3ex}\right],\end{eqnarray} \tag{ A.10 }$$

where ${\bar{n}}_{\pm }={\bar{n}}_{1}\pm {\bar{n}}_{2}.$ For the relevant case of not too small values of the squeezing parameter r,E_N can be well approximated with its value at equal mean thermal phonon number ${\bar{n}}_{-}=0,$

$$\begin{eqnarray}&&{E}_{N}({t}_{p})\simeq 4r-\mathrm{ln}\left[{\bar{n}}_{+}+1\right],\end{eqnarray} \tag{ A.11 }$$

showing that at this interaction time, the entanglement between the MR can be very large, even if starting from a relatively hot state, by properly tuning the ratio G₂/G₁, i.e., the intensity of the two tones. This large mechanical entanglement is achieved when the condition ${{\rm{e}}}^{{\rm{i}}{\phi }_{p}}=-1$ is also satisfied for a given integer p. This is obtained for any odd p when Δ = 0, or by properly adjusting the value of ${\mathcal{G}}$ for a given ${\rm{\Delta }}\ne 0,$ i.e., if

$$\begin{eqnarray}&&{{\mathcal{G}}}_{p}^{2}={{\rm{\Delta }}}^{2}\displaystyle \frac{d(2p-d)}{4{(p-d)}^{2}}\;\;\;d\;{\rm{odd}},\;0\lt d\lt 2p,\;d\ne p.\end{eqnarray} \tag{ A.12 }$$

This dynamical scheme for the generation of CV mechanical entanglement is similar to the Bogoliubov scheme proposed in [45] for entangling two optical cavity modes. It is extremely hard however to use it for entangling two mechanical modes as in the present case, because the cavity decay rate is comparable to ${\mathcal{G}}$ and Δ in typical situations, thereby strongly affecting the ideal Hamiltonian dynamics described here.

The system dynamics is described by the following QLE

$$\begin{eqnarray}\dot{\hat{a}} & = & -[\kappa +{\rm{i}}({{\rm{\Delta }}}_{0}+{\omega }_{-})]\hat{a}-{\rm{i}}\left[{E}_{1}{{\rm{e}}}^{-{\rm{i}}{\omega }_{+}t}+{E}_{2}{{\rm{e}}}^{{\rm{i}}{\omega }_{+}t}\right]+\sqrt{2\kappa }{\hat{a}}^{{\rm{in}}}\\ & & -{\rm{i}}\left[{g}_{1}\left({\hat{b}}_{1}+{\hat{b}}_{1}^{\dagger }\right)+{g}_{2}\left({\hat{b}}_{2}+{\hat{b}}_{2}^{\dagger }\right)\right]\hat{a},\\ {\dot{\hat{b}}}_{j} & = & -(\displaystyle \frac{{\gamma }_{j}}{2}+{\rm{i}}{\omega }_{j}){\hat{b}}_{j}-{\rm{i}}{g}_{j}{\hat{a}}^{\dagger }\hat{a}+\sqrt{{\gamma }_{j}}{\hat{b}}_{j}^{{\rm{in}}}\ \;\;\;j=1,2,\end{eqnarray} \tag{ B.1 }$$

where, here, differently from the description used in section 2, we are representing the cavity field in a reference frame rotating at the frequency ω₀ − (ω₂ − ω₁)/2, and we have introduced the frequencies

$$\begin{eqnarray*}&&{\omega }_{\pm }=\displaystyle \frac{{\omega }_{2}\pm {\omega }_{1}}{2}.\end{eqnarray*}$$

The other parameters and operators are defined in the main text.

If we perform a time dependent displacement, for both cavity and mechanical degrees of freedom, of the form

$$\begin{eqnarray}\hat{a}(t) & = & \delta \hat{a}(t)+\alpha (t),\\ {\hat{b}}_{j}(t) & = & \delta {\hat{b}}_{j}(t)+{\beta }_{j}(t),\end{eqnarray} \tag{ B.2 }$$

the QLE reduce to the form

$$\begin{eqnarray}\dot{\delta \hat{a}} & = & -\{\kappa +{\rm{i}}({{\rm{\Delta }}}_{0}+{\omega }_{-})+2{\rm{i}}{g}_{1}{\rm{Re}}[{\beta }_{1}(t)]+2{\rm{i}}{g}_{2}{\rm{Re}}[{\beta }_{2}(t)]\}\delta \hat{a}\\ & & -{\rm{i}}A(t)+\sqrt{2\kappa }{\hat{a}}^{{\rm{in}}}\\ & & -{\rm{i}}\alpha (t)[{g}_{1}(\delta {\hat{b}}_{1}+\delta {\hat{b}}_{1}^{\dagger })+{g}_{2}(\delta {\hat{b}}_{2}+\delta {\hat{b}}_{2}^{\dagger })]\\ & & -{\rm{i}}\left[{g}_{1}\left(\delta {\hat{b}}_{1}+\delta {\hat{b}}_{1}^{\dagger }\right)+{g}_{2}\left(\delta {\hat{b}}_{2}+\delta {\hat{b}}_{2}^{\dagger }\right)\right]\delta \hat{a},\end{eqnarray} \tag{ B.3 }$$

$$\begin{eqnarray}{\dot{\delta \hat{b}}}_{j} & = & -(\displaystyle \frac{{\gamma }_{j}}{2}+{\rm{i}}{\omega }_{j})\delta {\hat{b}}_{j}-{\rm{i}}{B}_{j}(t)+\sqrt{{\gamma }_{j}}{\hat{b}}_{j}^{{\rm{in}}}\\ & & -{\rm{i}}{g}_{j}[\alpha (t)\delta {\hat{a}}^{\dagger }+\alpha {(t)}^{*}\delta \hat{a}]-{\rm{i}}{g}_{j}\delta {\hat{a}}^{\dagger }\delta \hat{a},\end{eqnarray} \tag{ B.4 }$$

where the new driving terms, A(t) and B_j(t) read

$$\begin{eqnarray}A(t) & = & -{\rm{i}}\left[{E}_{1}{{\rm{e}}}^{-{\rm{i}}{\omega }_{+}t}+{E}_{2}{{\rm{e}}}^{{\rm{i}}{\omega }_{+}t}\right]-\displaystyle \frac{\partial \alpha (t)}{\partial t}-[\kappa +{\rm{i}}({{\rm{\Delta }}}_{0}+{\omega }_{-})]\alpha (t)\\ & & -2{\rm{i}}({g}_{1}{\rm{Re}}[{\beta }_{1}(t)]+{g}_{2}{\rm{Re}}[{\beta }_{2}(t)])\alpha (t),\\ {B}_{j}(t) & = & -\displaystyle \frac{\partial {\beta }_{j}(t)}{\partial t}-(\displaystyle \frac{{\gamma }_{j}}{2}+{\rm{i}}{\omega }_{j}){\beta }_{j}(t)-{\rm{i}}{g}_{j}{| \alpha (t)| }^{2}.\end{eqnarray} \tag{ B.5 }$$

When g₁ and g₂ are sufficiently small and we chose α(t) and β_j(t) such that A(t) = 0 and B_j(t) = 0, then the nonlinear terms, i.e. the last terms in the two equations (B.3) and (B.4), can be neglected. The equations A(t) = 0 and B_j(t) = 0 define a set of nonlinear differential equations with periodic driving for the parameters α(t) and β_j(t). The solution can be evaluated perturbatively in the small parameters g₁ and g₂[46]. Here we assume g₁ = g₂ ≡ g and we observe that the solutions for α(t) and β_j(t), with initial condition α(0) = β_j(0) = 0, contain, respectively, only even and odd powers of g,

$$\begin{eqnarray}\alpha (t) & = & \displaystyle \sum _{\begin{array}{c}p=0\\ p\ {\rm{even}}\end{array}}^{\infty }{g}^{p}{\alpha }^{(p)}(t),\\ {\beta }_{j}(t) & = & \displaystyle \sum _{\begin{array}{c}p=1\\ p\ {\rm{odd}}\end{array}}^{\infty }{g}^{p}{\beta }_{j}^{(p)}(t).\end{eqnarray} \tag{ B.6 }$$

The equations for each component of these expansions can be written in the form

$$\begin{eqnarray}{\dot{\alpha }}^{(p)}(t) & = & -z\;{\alpha }^{(p)}(t)+{{\rm{\Xi }}}_{\alpha }^{(p)}(t),\\ {\dot{\beta }}_{j}^{(p)}(t) & = & -{w}_{j}\;{\beta }_{j}^{(p)}(t)+{{\rm{\Xi }}}_{\beta }^{(p)}(t),\end{eqnarray} \tag{ B.7 }$$

where

$$\begin{eqnarray}z & = & \kappa +{\rm{i}}({{\rm{\Delta }}}_{0}+{\omega }_{-}),\\ {w}_{j} & = & \displaystyle \frac{{\gamma }_{j}}{2}+{\rm{i}}{\omega }_{j},\end{eqnarray} \tag{ B.8 }$$

and the driving terms are defined recursively as

$$\begin{eqnarray}{{\rm{\Xi }}}_{\alpha }^{(p)}(t) & = & -2{\rm{i}}\displaystyle \sum _{q=0}^{p-1}{\alpha }^{(q)}(t)\\ & & \times ({\rm{Re}}[{\beta }_{1}^{(p-q-1)}(t)]+{\rm{Re}}[{\beta }_{2}^{(p-q-1)}(t)]),\\ {{\rm{\Xi }}}_{\beta }^{(p)}(t) & = & -{\rm{i}}\displaystyle \sum _{q=0}^{p-1}{\alpha }^{(q)}(t)\ {\alpha }^{(p-q-1)}{(t)}^{*},\end{eqnarray} \tag{ B.9 }$$

with the initial condition

$$\begin{eqnarray}{{\rm{\Xi }}}_{\alpha }^{(0)}(t) & = & {E}_{1}\;{{\rm{e}}}^{-{\rm{i}}\;{\omega }_{+}\;t}+{E}_{2}\;{{\rm{e}}}^{{\rm{i}}\;{\omega }_{+}\;t},\\ {{\rm{\Xi }}}_{\beta }^{(0)}(t) & = & 0.\end{eqnarray} \tag{ B.10 }$$

In particular they can always be rewritten as sums of exponential functions of the form

$$\begin{eqnarray}{{\rm{\Xi }}}_{\alpha }^{(p)}(t) & = & \displaystyle \sum _{n}{\chi }_{\alpha }^{(p,n)}\;{{\rm{e}}}^{{\zeta }_{\alpha }^{(p,n)}\;t},\\ {{\rm{\Xi }}}_{\beta }^{(p)}(t) & = & \displaystyle \sum _{n}{\chi }_{\beta }^{(p,n)}\;{{\rm{e}}}^{{\zeta }_{\beta }^{(p,n)}\;t},\end{eqnarray} \tag{ B.11 }$$

with ${\chi }_{\alpha }^{(p,n)}{\chi }_{\beta }^{(p,n)}{\zeta }_{\alpha }^{(p,n)}$ and ${\zeta }_{\beta }^{(p,n)}$ time-independent complex coefficients, whose specific form can be computed iteratively. Moreover, the expression for ${\alpha }^{(p)}(t)$ and ${\beta }_{j}^{(p)}(t)$ are found integrating equation (B.7) and are given by

$$\begin{eqnarray}{\alpha }^{(p)}(t) & = & \displaystyle \sum _{n}\displaystyle \frac{{\chi }_{\alpha }^{(p,n)}}{z+{\zeta }_{\alpha }^{(p,n)}}({{\rm{e}}}^{{\zeta }_{\alpha }^{(p,n)}t}-{{\rm{e}}}^{-zt}),\\ {\beta }_{j}^{(p)}(t) & = & \displaystyle \sum _{n}\displaystyle \frac{{\chi }_{\beta }^{(p,n)}}{{w}_{j}+{\zeta }_{\beta }^{(p,n)}}({{\rm{e}}}^{{\zeta }_{\beta }^{(p,n)}t}-{{\rm{e}}}^{-{w}_{j}t}).\end{eqnarray} \tag{ B.12 }$$

We note that all the coefficients ${\zeta }_{\alpha }^{(p,n)}$ and ${\zeta }_{\beta }^{(p,n)}$ have non-positive real parts, ${\rm{Re}}[{\zeta }_{\alpha }^{(p,n)}],{\rm{Re}}[{\zeta }_{\beta }^{(p,n)}]\leqslant 0,$ thus the large-time solutions ${\alpha }^{(p)}(t\to \infty )\equiv {\alpha }_{\mathrm{st}}^{(p)}(t)$ and ${\beta }_{j}^{(p)}(t\to \infty )\equiv {\beta }_{j,\mathrm{st}}^{(p)}(t)$ are found from equation (B.12) by keeping only the terms for which ${\zeta }_{\alpha }^{(p,n)}$ and ${\zeta }_{\beta }^{(p,n)}$ are purely imaginary, that can be shown to be equal to ${\rm{i}}\;n\;{\omega }_{+}$ with n odd and even integer respectively. In particular ${\alpha }_{\mathrm{st}}^{(p)}(t)$ and ${\tilde{\beta }}_{j,\mathrm{st}}^{(p)}(t)$ are periodic functions (with period 2π/ω₊) which contains frequency components that are, respectively, odd and even multiples of ω₊,

$$\begin{eqnarray}{\alpha }_{\mathrm{st}}^{(p)}(t) & = & \displaystyle \sum _{\begin{array}{c}n=-p-1\\ n\ {\rm{odd}}\end{array}}^{p+1}\displaystyle \frac{{\tilde{\chi }}_{\alpha }^{(p,n)}}{z+{\rm{i}}\;n\;{\omega }_{+}}{{\rm{e}}}^{{\rm{i}}n{\omega }_{+}t},\\ {\beta }_{j,\mathrm{st}}^{(p)}(t) & = & \displaystyle \sum _{\begin{array}{c}n=-p-1\\ n\ {\rm{even}}\end{array}}^{p+1}\displaystyle \frac{{\tilde{\chi }}_{\beta }^{(p,n)}}{{w}_{j}+{\rm{i}}\;n\;{\omega }_{+}}{{\rm{e}}}^{{\rm{i}}n{\omega }_{+}t},\end{eqnarray} \tag{ B.13 }$$

where z and w_j are defined in equation (B.8), and ${\tilde{\chi }}_{\alpha }^{(p,n)}(t)$ and ${\tilde{\chi }}_{\beta }^{(p,n)}(t)$ are the coefficients that correspond to those particular parameters ${\zeta }_{\alpha }^{(p,n)}$ and ${\zeta }_{\beta }^{(p,n)}$ that are imaginary.

B.1. Resonant and non-resonant terms

The QLE, in the interaction picture with respect to the Hamiltonian ${\hat{H}}_{0}={\hslash }({\omega }_{-}{a}^{\dagger }a+{\omega }_{1}{\hat{b}}_{1}^{\dagger }{\hat{b}}_{1}+{\omega }_{2}{\hat{b}}_{2}^{\dagger }{\hat{b}}_{2}),$ reduce to

$$\begin{eqnarray}\dot{\delta \hat{a}} & = & -(\kappa +{\rm{i}}{{\rm{\Delta }}}_{0})\delta \hat{a}+\sqrt{2\kappa }{\hat{a}}^{{\rm{in}}}\\ & & +2{\rm{i}}({g}_{1}{\rm{Re}}[{\beta }_{1}(t)]+{g}_{2}{\rm{Re}}[{\beta }_{2}(t)])\delta \hat{a}\\ & & -{\rm{i}}\alpha (t){{\rm{e}}}^{{\rm{i}}{\omega }_{-}t}\left[{g}_{1}(\delta {\hat{b}}_{1}{{\rm{e}}}^{-{\rm{i}}{\omega }_{1}t}+\delta {\hat{b}}_{1}^{\dagger }{{\rm{e}}}^{{\rm{i}}{\omega }_{1}t})\right.\\ & & \left.+{g}_{2}(\delta {\hat{b}}_{2}{{\rm{e}}}^{-{\rm{i}}{\omega }_{2}t}+\delta {\hat{b}}_{2}^{\dagger }{{\rm{e}}}^{{\rm{i}}{\omega }_{2}t})\right],\\ {\dot{\delta \hat{b}}}_{j} & = & -\displaystyle \frac{{\gamma }_{j}}{2}\delta {\hat{b}}_{j}+\sqrt{{\gamma }_{j}}{\hat{b}}_{j}^{{\rm{in}}}-{\rm{i}}{g}_{j}{{\rm{e}}}^{{\rm{i}}{\omega }_{j}t}[{{\rm{e}}}^{{\rm{i}}{\omega }_{-}t}\alpha (t)\delta {\hat{a}}^{\dagger }+{\rm{h.c.}}].\end{eqnarray} \tag{ B.14 }$$

Before proceeding, we note that we can include the dc component of β_j(t) into the cavity detuning, hence we introduce ${\rm{\Delta }}={{\rm{\Delta }}}_{0}+2{\displaystyle \sum }_{j=\mathrm{1,2}}{g}_{j}{\rm{Re}}[{\sum }_{p}{g}^{p}\frac{{\tilde{\chi }}_{\beta }^{(p,0)}}{{w}_{j}}],$ according to the notation introduced in equation (8),

$$\begin{eqnarray}&&\displaystyle \sum _{p}{g}^{p}\displaystyle \frac{{\tilde{\chi }}_{\beta }^{(p,0)}}{{w}_{j}}\equiv {\beta }_{j}^{{\rm{dc}}}.\end{eqnarray} \tag{ B.15 }$$

Moreover we can isolate the resonant terms of the QLE, namely the terms with time-independent coefficients, by considering the lowest order frequency components of α(t), i.e.

$$\begin{eqnarray}&&{\alpha }_{\pm }=\displaystyle \sum _{p}{g}^{p}\displaystyle \frac{{\tilde{\chi }}_{\alpha }^{(p,\pm 1)}}{z\pm {\rm{i}}\;{\omega }_{+}},\end{eqnarray} \tag{ B.16 }$$

corresponding to the frequencies ±ω₊, and defining

$$\begin{eqnarray}{\bar{\beta }}_{j}(t) & = & {\beta }_{j}(t)-{\beta }_{j}^{{\rm{dc}}},\\ {\bar{\alpha }}_{+}(t) & = & \alpha (t)-{{\rm{e}}}^{{\rm{i}}{\omega }_{+}t}{\alpha }_{+},\end{eqnarray} \tag{ B.17 }$$

$$\begin{eqnarray}&&{\bar{\alpha }}_{-}(t)=\alpha (t)-{{\rm{e}}}^{-{\rm{i}}{\omega }_{+}t}{\alpha }_{-}.\end{eqnarray} \tag{ B.18 }$$

Thereby we find

$$\begin{eqnarray}\dot{\delta \hat{a}} & = & -(\kappa +{\rm{i}}{\rm{\Delta }})\delta \hat{a}-{\rm{i}}{\alpha }_{-}{g}_{1}\delta {\hat{b}}_{1}^{\dagger }-{\rm{i}}{\alpha }_{+}{g}_{2}\delta {\hat{b}}_{2}\\ & & +\sqrt{2\kappa }{\hat{a}}^{{\rm{in}}}+{F}_{a}(t),\\ {\dot{\delta \hat{b}}}_{1} & = & -\displaystyle \frac{{\gamma }_{1}}{2}\delta {\hat{b}}_{1}-{\rm{i}}{g}_{1}{\alpha }_{-}\delta {\hat{a}}^{\dagger }+\sqrt{{\gamma }_{1}}{\hat{b}}_{1}^{{\rm{in}}}+{F}_{{b}_{1}}(t),\\ {\dot{\delta \hat{b}}}_{2} & = & -\displaystyle \frac{{\gamma }_{2}}{2}\delta {\hat{b}}_{2}-{\rm{i}}{g}_{2}{{\alpha }_{+}}^{*}\delta \hat{a}+\sqrt{{\gamma }_{2}}{\hat{b}}_{2}^{{\rm{in}}}+{F}_{{b}_{2}}(t),\end{eqnarray} \tag{ B.19 }$$

where ${F}_{a}(t),{F}_{{b}_{1}}(t)$ and ${F}_{{b}_{2}}(t)$ account for the terms with time-dependent coefficients and are given by

$$\begin{eqnarray}{F}_{a}(t) & = & +2{\rm{i}}({g}_{1}{\rm{Re}}[{\bar{\beta }}_{1}(t)]+{g}_{2}{\rm{Re}}[{\bar{\beta }}_{2}(t)])\delta \hat{a}\\ & & -{\rm{i}}\alpha (t){{\rm{e}}}^{{\rm{i}}{\omega }_{-}t}({g}_{1}\delta {\hat{b}}_{1}{{\rm{e}}}^{-{\rm{i}}{\omega }_{1}t}+{g}_{2}\delta {\hat{b}}_{2}^{\dagger }{{\rm{e}}}^{{\rm{i}}{\omega }_{2}t})\\ & & -{\rm{i}}{\bar{\alpha }}_{-}(t){{\rm{e}}}^{{\rm{i}}{\omega }_{-}t}{g}_{1}\delta {\hat{b}}_{1}^{\dagger }{{\rm{e}}}^{{\rm{i}}{\omega }_{1}t}-{\rm{i}}{\bar{\alpha }}_{+}(t){{\rm{e}}}^{{\rm{i}}{\omega }_{-}t}{g}_{2}\delta {\hat{b}}_{2}{{\rm{e}}}^{-{\rm{i}}{\omega }_{2}t},\\ {F}_{{b}_{1}}(t) & = & -{\rm{i}}{g}_{1}{{\rm{e}}}^{{\rm{i}}{\omega }_{1}t}[{{\rm{e}}}^{{\rm{i}}{\omega }_{-}t}{\bar{\alpha }}_{-}(t)\delta {\hat{a}}^{\dagger }+{{\rm{e}}}^{-{\rm{i}}{\omega }_{-}t}{\alpha }^{*}(t)\delta \hat{a}],\\ {F}_{{b}_{2}}(t) & = & -{\rm{i}}{g}_{2}{{\rm{e}}}^{{\rm{i}}{\omega }_{2}t}[{{\rm{e}}}^{{\rm{i}}{\omega }_{-}t}\alpha (t)\delta {\hat{a}}^{\dagger }+{{\rm{e}}}^{-{\rm{i}}{\omega }_{-}t}{\bar{\alpha }}_{+}{(t)}^{*}\delta \hat{a}].\end{eqnarray} \tag{ B.20 }$$

In particular we can introduce the linearized coupling strength ${G}_{1}^{(\mathrm{tot})}=g{\alpha }_{-}$ and ${G}_{2}^{(\mathrm{tot})}=g{\alpha }_{+},$ with α_± defined in equation (B.16). The expressions introduced in equation (7) correspond to the expansion of these parameters at zeroth order in g (see also equation (B.22)).

We are interested in the regime in which the terms in equations (B.20) with time-dependent coefficients are negligible. They can be neglected when ${g}_{j}\;| {\alpha }_{\mathrm{st}}(t)| ,{g}_{j}\;| {\bar{\beta }}_{j,\mathrm{st}}(t)| ,\kappa \ll {\rm{min}}\{{\omega }_{1},{\omega }_{2},| {\omega }_{1}-{\omega }_{2}| \}.$ In particular this condition is true when it is valid for the lowest order term in the expansion in power of g. In details, the non resonant terms can be neglected when

$$\begin{eqnarray}&&g\;| {\alpha }_{\pm }| ,\kappa \ll {\rm{min}}\{{\omega }_{1},{\omega }_{2},| {\omega }_{1}-{\omega }_{2}| \}.\end{eqnarray} \tag{ B.21 }$$

When this condition is fulfilled the parameters α(t) and β_j(t) can be safely expanded at the lowest order in g. Specifically they can be approximated as

$$\begin{eqnarray}{\alpha }_{-} & \simeq & \displaystyle \frac{-{\rm{i}}{E}_{1}}{z-{\rm{i}}{\omega }_{+}},\ \ \ \ \ \ \ \ \ {\alpha }_{+}\simeq \displaystyle \frac{-{\rm{i}}{E}_{2}}{z+{\rm{i}}{\omega }_{+}},\\ {\alpha }_{\mathrm{st}}(t) & \simeq & {\alpha }_{\mathrm{st}}^{(0)}(t)={\alpha }_{-}{{\rm{e}}}^{-{\rm{i}}{\omega }_{+}t}+{\alpha }_{+}{{\rm{e}}}^{{\rm{i}}{\omega }_{+}t},\\ {\beta }_{j}^{{\rm{dc}}} & \simeq & \displaystyle \frac{-{\rm{i}}{g}_{j}}{{w}_{j}}[{\alpha }_{-}{\alpha }_{-}^{*}+{\alpha }_{+}{\alpha }_{+}^{*}],\\ {\bar{\beta }}_{j,\mathrm{st}}(t) & \simeq & {g}_{j}\;{\beta }_{j,\mathrm{st}}^{(1)}(t)-{\beta }_{j}^{{\rm{dc}}}\\ & = & -{\rm{i}}{g}_{j}[\displaystyle \frac{{\alpha }_{-}{\alpha }_{+}^{*}}{{w}_{j}-2{\rm{i}}{\omega }_{+}}{{\rm{e}}}^{-2{\rm{i}}\;{\omega }_{+}\;t}\displaystyle \frac{{\alpha }_{+}{\alpha }_{-}^{*}}{{w}_{j}+2{\rm{i}}{\omega }_{+}}{{\rm{e}}}^{2{\rm{i}}\;{\omega }_{+}\;t}].\end{eqnarray} \tag{ B.22 }$$

Moreover the parameters ${\bar{\alpha }}_{\pm }(t)$ defined in equation (B.18) are zero. Using these expressions the QLE in equation (B.19) can be rewritten as

$$\begin{eqnarray}\dot{\delta \hat{a}} & = & -(\kappa +{\rm{i}}{\rm{\Delta }})\delta \hat{a}-{\rm{i}}{G}_{1}\delta {\hat{b}}_{1}^{\dagger }-{\rm{i}}{G}_{2}\delta {\hat{b}}_{2}\\ & & +\sqrt{2\kappa }{\hat{a}}^{{\rm{in}}}+{F}_{a}(t),\\ {\dot{\delta \hat{b}}}_{1} & = & -\displaystyle \frac{{\gamma }_{1}}{2}\delta {\hat{b}}_{1}-{\rm{i}}{G}_{1}\delta {\hat{a}}^{\dagger }+\sqrt{{\gamma }_{1}}{\hat{b}}_{1}^{{\rm{in}}}+{F}_{{b}_{1}}(t),\\ {\dot{\delta \hat{b}}}_{2} & = & -\displaystyle \frac{{\gamma }_{2}}{2}\delta {\hat{b}}_{2}-{\rm{i}}{G}_{2}^{*}\delta \hat{a}+\sqrt{{\gamma }_{2}}{\hat{b}}_{2}^{{\rm{in}}}+{F}_{{b}_{2}}(t),\end{eqnarray} \tag{ B.23 }$$

with G₁ and G₂ defined in equation (7) and

$$\begin{eqnarray*}{F}_{a}(t) & \simeq & +2{\rm{i}}({g}_{1}^{2}{\rm{Re}}[{\bar{\beta }}_{1,\mathrm{st}}^{(1)}(t)]+{g}_{2}^{2}{\rm{Re}}[{\bar{\beta }}_{2,\mathrm{st}}^{(1)}(t)])\delta \hat{a}\\ & & -{\mathrm{ie}}^{{\rm{i}}\;{\omega }_{-}\;t}{\alpha }_{\mathrm{st}}^{(0)}(t)({g}_{1}\delta {\hat{b}}_{1}{{\rm{e}}}^{-{\rm{i}}{\omega }_{1}t}+{g}_{2}\delta {\hat{b}}_{2}^{\dagger }{{\rm{e}}}^{{\rm{i}}{\omega }_{2}t}),\end{eqnarray*}$$

$$\begin{eqnarray}{F}_{{b}_{1}}(t) & \simeq & -{\rm{i}}{g}_{1}{{\rm{e}}}^{{\rm{i}}}({\omega }_{1}+{\omega }_{-})t{{\alpha }_{\mathrm{st}}^{(0)}}^{*}(t)\delta \hat{a},\\ {F}_{{b}_{2}}(t) & \simeq & -{\rm{i}}{g}_{2}{{\rm{e}}}^{{\rm{i}}}({\omega }_{2}+{\omega }_{-})t{\alpha }_{\mathrm{st}}^{(0)}(t)\delta {\hat{a}}^{\dagger }.\end{eqnarray} \tag{ B.24 }$$

When the time-dependent coefficients are neglected these equations reduce to equations (4)–(6).