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

We investigate a general scheme for generating, either dynamically or in the steady state, continuous variable entanglement between two mechanical resonators with different frequencies. We employ an optomechanical system in which a single optical cavity mode driven by a suitably chosen two-tone field is coupled to the two resonators. Significantly large mechanical entanglement can be achieved, which is extremely robust with respect to temperature.


I. INTRODUCTION
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 Ref. [2], while entanglement between two Josephson-junction qubits has been detected in Refs. [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 Ref. [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 singlephonon 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][12][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. Refs. [16][17][18] considered the steady state of different systems of driven cavities: Ref. [16] focused on two mirrors of a ring cavity, while Ref. [17] assumed to drive two independent linear cavities with two-mode squeezed light transferring its entanglement to the cavity end-mirrors. Ref. [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][23][24][25][26][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][30][31][32][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][35][36][37][38][39][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 Refs. [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 II we derive the effective quantum Langevin equations (QLE) describing the dynamics of the system in the RWA. In Section III 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 IV 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 V 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 VI 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.

II. SYSTEM HAMILTONIAN AND DERIVATION OF THE EFFECTIVE LANGEVIN EQUATIONS
As shown in Fig. 1, we consider an optical cavity mode with resonance frequency ω c and annihilation operatorâ interacting via the usual optomechanical interaction with two different MRs, with frequencies ω 1 and ω 2 and annihilation operatorsb 1 andb 2 respectively. The cavity mode is bichromatically driven at the two frequencies ω 0 + ω 1 and ω 0 − ω 2 , with the reference frequency ω 0 detuned from the cavity resonance by a quantity ∆ 0 = ω c − ω 0 . If we describe the cavity field in a reference frame rotating at the frequency ω 0 , then the system Hamiltonian is given bŷ This means that the cavity mode is simultaneously driven on the blue sideband associated with the MR with annihilation operatorb 1 , and on the red sideband associated with the MR withb 2 . The nonzero detuning ∆ 0 makes the present scheme different from the one studied in the supplementary material of Ref. [34] which restricts to the resonant case ∆ 0 = 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.
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 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 δâ(t) and δb j (t) can be derived, in the interaction picture with respect to the HamiltonianĤ 0 = ω 1b † 1b 1 + ω 2b † 2b 2 , by neglecting the non-resonant/timedependent components of the effective linearized interactions. It is possible to prove that this approach is justified when (see Eq. (B21)) The corresponding QLE including thermal noise and dissipation at rates κ and γ j for the cavity and the mechanical mode j ∈ 1, 2 respectively, are where are the (generally complex) linear optomechanical couplings, andâ in andb in j are standard input noise operators with zero mean, whose only nonzero correlation functions is the mean thermal phonon number of the j-th 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 β DC j (see The cavity mode is bichromatically driven at the two frequencies ω 0 + ω 1 and ω 0 − ω 2 . Large and robust entanglement of the two mechanical resonators can be generated either dynamically or in the steady state. Two weak probe fields with detuning ∆ p j = ω c j − ω p j = ω j , j = 1, 2, are then sent into the cavity. By homodyning the probe mode outputs, the mechanical quadratures (x j , p j ) are therefore measured, which allows one to construct the correlation matrix of the quadratures from which entanglement can be derived in a straightforward way.

Appendix B). Specifically
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 Eqs. (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, γ 1 = γ 2 = γ. In such a case, the system is stable if and only if This stability condition reduces to the one derived in the supplementary material of Ref. [34] in the case ∆ = 0. We see that a nonzero detuning generally helps in keeping the system stable.

III. DARK AND BRIGHT BOGOLIUBOV MODES
The coherent dynamics corresponding to the Eqs. (4)-(6), is described by the effective linearized Hamiltonian 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 1 and G 2 real. Eq. (10) naturally suggests to introduce two effective mechanical modes allowing to simplify the system dynamics. We assume for the moment G 2 > G 1 , which is a sufficient condition for stability (see Eq. (9)), and definê where Eqs. (11)-(12) define a Bogoliubov unitary transformation of the mechanical mode operators, which can also be written aŝ withŜ (r) the two-mode squeezing operator. The Bogoliubov modeβ 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β 2 is the "bright" mode interacting with the cavity mode. This is equivalent to say that the dark modeβ 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β 2 and δâ. The Bogoliubov mode description has been already employed in cavity optomechanics, associated to two optical modes in Refs. [34,45], and to two mechanical modes in Refs. [36,38] (see Appendix A for a derivation of the normal modes of the system and a study of its Hamiltonian dynamics).

A. 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 γ 1 = γ 2 ≡ γ assume the simple form whereβ in j , j = 1, 2, are two correlated thermal noise operators whose only nonzero correlation functions are with the effective mean thermal phonon numbers n eff 1 (r) =n 1 cosh 2 r + (n 2 + 1) sinh 2 r, n eff 2 (r) =n 2 cosh 2 r + (n 1 + 1) sinh 2 r, (20) and the inter-mode correlation m(r) = cosh r sinh r (n 1 +n 2 + 1) .
If γ 1 γ 2 a dissipative coupling term between the two Bogoliubov modes appears, which however does not have relevant effects because it is proportional to |γ 1 − γ 2 | which is typically very small with respect to all other damping rates. The dynamics associated with Eqs. (15)-(17) is simple: the bright mechanical modeβ 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 β = β 1 ,β 2 ,β 1 † ,β 2 † , whose elements are C β j,k = {β} j {β} k is given, at the steady state, by with the number of excitation of the cooled bright mode and the correlations between the two Bogoliubov modes respectively given bȳ where 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 1 and b 2 by inverting the Bogoliubov transformation introduced in Eqs. (11) and (12). The result is and where now n b1 =n eff 1 + sinh 2 r 1 +n eff 1 +n cool 2 −2 cosh r sinh r Re m β , n b2 =n cool 2 + sinh 2 r 1 +n eff 1 +n cool 2 −2 cosh r sinh r Re m β , The entanglement between modes b 1 and b 2 , measured by means of the logarithmic negativity [47,48], can be easily expressed in terms of these matrix elements as [49] When the collective cooperativity C − is sufficiently large, i.e., C − m(r), thenm β (r) is negligible (see Eq. (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 (n cool 2 (r) n eff 2 (r)), corresponding to an entangled state for the original mechanical modes. In this case the steady state correlation matrix for the Bogoliubov modes, in Eq. (22), reduces to the correlation matrix of a state given by the product of two thermal states with occupanciesn eff 1 (r) andn cool 2 (r) respectively. For the two MR of interest, associated with the operatorb 1 andb 2 , such a state is just a two-mode squeezed thermal state [50] whereŜ (r) is given in Eq. (14), and is the density matrix of the thermal equilibrium state of a resonator with occupancyn. 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 Fig. 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 Eqs. (4)- (6). Figure 2 refers to an experimentally achievable set of parameters, γ = 10 s −1 , κ = 10 5 s −1 , G 2 = 10 5 s −1 , ∆ = 10 3 s −1 , and to different values of mean thermal phonon numbers n 1 ,n 2 , and of the ratio G 1 /G 2 . 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 forn 1 = 2000,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 (κ 2 + ∆ 2 )/(G 2 κ) (see Eqs. (15)-(17)).
Eq. (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 − 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 1 and G 2 , is far from being trivial. In fact, r increases when G 1 → G 2 , which however implies, at a fixed value of G 2 , a decreasing value of G and therefore of C − (see Eq. (13) and Eq. (25)); moreover increasing r has also the unwanted effect of increasingm β (r) that is the correlations between the two Bogoliubov modes (see Eqs. (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 occupanciesn 1 andn 2 . At a given value of G 1 , this is obtained by taking a sufficiently large value of the associated singlemode cooperativity, C 1 = 2G 2 1 /κγ 1, and correspondingly optimizing the value of G 2 , i.e., of r. In fact, the logarithmic negativity associated with the stationary state of Eq. (32) can be evaluated in terms of the parameter − n − (r) 2 + 4 [n + (r) + 1] 2 sinh 2 r cosh 2 r , wheren ± (r) =n eff 1 (r) ±n cool 2 (r), andn cool 2 (r) can be explicitly rewritten in terms of the cooperativity C 1 as n cool 2 (r) = n 2 cosh 2 r + (n 1 + 1) sinh 2 r The dependence of E N versus r, for given values of C 1 ,n 1 and n 2 , shows a maximum and then decays to zero for large r (see Fig. 3 which refers to C 1 = 2 × 10 4 andn 1 = 200,n 2 = 100). This behavior is described by a very simple approximated expression valid in the limit C 1 e 2r e −2r , with not very largen 1,2 , and when δ, → 0 (corresponding to γ, ∆ κ), which exhibits a minimum (hence corresponding to maximum entanglement) as a function of r at given by ν opt ∼ 2(1+n 1 +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 r opt , G is reduced and the cooling dynamics becomes slow as compared to the standard mechanical dissipation, which takes place at rate ∼ γ (n j + 1), so that the correlations between the MRs cannot be efficiently generated. On the other hand, at small r r 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. Fig. 3 also shows that the simplified expression of Eq. (36) provides a simple but valid approximation for large C 1 and a very good estimate of the optimal value of the twomode squeezing parameter r, i.e., of G 1 /G 2 , given by Eq. (37). The corresponding value of the logarithmic negativity is and shows that once that the ratio G 1 /G 2 is optimized, the achievable stationary entanglement between the two MRs increases with increasing C 1 /(n 1 +n 2 ). The above analysis of the stationary entanglement of the two MRs extends the results of Ref. [34] in various directions. First of all, our model extends to the case of nonzero detuning ∆ a model discussed in the Supplementary material of Ref. [34]. We see that a nonzero detuning has a limited effect of the dynamic of entanglement generation, providing only an effective increase ofn cool 2 , which however becomes negligible as soon as ∆ κ (see Eq. (35)). Moreover, Ref. [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,n 1 andn 2 .

IV. DYNAMICAL EVOLUTION IN THE CASE OF EQUAL COUPLINGS
In the special case of equal couplings G 1 = G 2 ≡ G, i.e., 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 Ref. [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 Eqs. (4)-(6) are the same as those studied in Ref. [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 Eq. (10) in terms of mechanical and optical quadratures, using the expressions δb j = (x j + ip j )/ √ 2, j = 1, 2, and δâ = (X + iŶ)/ wherex ± = (x 1 ±x 2 )/ √ 2,p ± = (p 1 ±p 2 )/ √ 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 thatx + andp − are two commuting conserved observables. One gets (see also Refs. [43,44,51]) 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, . . ., wherê 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 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 ∆ κ. 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 Fig. 4, which refers to the parameter set γ = 10 s −1 , κ = 10 5 s −1 , G = 10 5 s −1 , n 1 = 200,n 2 = 100, and to three different values of the detuning, ∆ = 10 3 s −1 (black dashed line), ∆ = 10 4 s −1 (red full line), and ∆ = 10 5 s −1 (blue full line). We see that an appreciable value of E N (even though smaller than the one achievable at the samen 1 andn 2 after the optimization of G 1 /G 2 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.

V. EFFECT OF THE COUNTER-ROTATING TERMS. STUDY OF THE EXACT DYNAMICS
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 Ref. [34]). It is therefore interesting to study their effect by comparing the above predictions, both in the case of G 2 > G 1 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 Fig. 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 Eq. (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.
In Fig. 5 we compare the time evolution of the entanglement evaluated with and without the time-dependent terms when G 2 > G 1 . The parameters used in these plots are consistent with those used in Fig. 2. Specifically the three red curves in Figs. 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 Fig. 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 Eq. (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  (B23), which takes into account the non-resonant terms by considering the expansion of the steady state solutions of α(t) and β j (t), at first order in g as defined in Eq. (B22). The blue lines are evaluated instead with Eq. (B19) by considering the expansion for α(t) and β j (t), calculated iteratively with Eqs. (B6), (B9)-(B12), up to sixth order in g; in particular these results take into account the full dynamics of the average fields α(t) and β j (t), with initial condition α(0) = β j (0) = 0, and not only the steady state as in the case of the green lines. The solid lines refer to ω 2 = 100κ and ω 1 = 50κ, while the dashed lines refer to ω 2 = 50κ and ω 1 = 25κ. The other parameters are G 2 = κ, ∆ = 0.01κ, γ = 10 −4 κ, and g = 10 −4 κ.  Fig. 5. The solid lines are found for ∆ = 0.01κ and the dashed lines for ∆ = 5κ. In (a) γ = 0.03κ, ω 1 = 58κ, ω 2 = 100κ,n 1 = 2 andn 2 = 1; in (b) γ = 0.01κ, ω 1 = 58κ, ω 2 = 100κ,n 1 = 2 and n 2 = 1; in (c) γ = 0.001κ, ω 1 = 51κ, ω 2 = 100κ,n 1 = 20 andn 2 = 10. The other parameters are G 1 = G 2 = κ and g = 10 −4 κ.
Eq. (3), in order to eliminate the fast rotating terms, the ratios ω j /G j have to be much larger than one. Although the dashedblue curves are evaluated for a ratio ω 1 /G 1 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 Fig. 6 we study the case of equal couplings G 1 = G 2 . 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 Fig. 6 (a). Larger discrepancies are found when the mechanical dissipation is reduced as in Fig. 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 behaviour 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 Eq. (9)). We have confirmed the unstable nature of the timedependent 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 Fig. 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 Fig. 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 (ω j 10 2 κ) (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.

VI. STRATEGIES FOR THE EXPERIMENTAL DETECTION OF MECHANICAL ENTANGLEMENT
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 Ref. [53] and then implemented in the electromechanical entanglement experiment of Ref. [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 ∆ p j = ω j , j = 1, 2, the probe mode adiabatically follows the MR dynamics, and the output of the readout cavity a out j is given by (see Fig. 1) [53] a out with G p j 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 .

A. Concluding remarks
We have studied in detail a general scheme for the generation of large and robust CV entanglement between two MRs with different frequencies through their coupling with a single, bichromatically driven cavity mode. The scheme extends and generalizes in various directions similar schemes exploiting driven cavity modes [34,36,38,44] for entangling two MRs or two cavity modes. The scheme is able to generate a remarkably large entanglement between two macroscopic oscillators in the stationary state, i.e., with virtually infinite lifetime, and it is quite robust because one can achieve appreciably large CV entanglement even with thermal occupancies of the order of 10 3 . The scheme is particularly efficient in the limit where counter-rotating terms due to the bichromatic driving of the cavity mode are negligible, and we have verified with a careful numerical analysis that this is well justified when the two mechanical frequencies are sufficiently large ω j 10 2 κ.

VII. ACKNOWLEDGMENTS
This work has been supported by the European Commission (ITN-Marie Curie project cQOM and FET-Open Project iQUOEMS), by MIUR (PRIN 2011). j = 1, 2, and β 1 (t) = β 1 (0), one gets 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 Eqs. (A4)-(A5) show that it occurs when sin∆t/2 = 0, i.e., t p = 2πp/ ∆ 2 + 4G 2 , p = 1, 2, . . .. At these time instants one has δb 1 (t p ) = cosh 2 r − e iφ p sinh 2 r δb 1 (0) (A6) where φ p = πp 1 + ∆/∆ . In particular, if e iφ p = −1 one gets δb † 2 (t p ) = − cosh 2rδb † 2 (0) − sinh 2rδb 1 (0), (A9) 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,Ŝ (2r) (see Eq. (14)) to their initial state. In the usual case of an initial thermal state for the two MRs with mean thermal phonon numbersn j , the state at time t p is therefore a two-mode squeezed thermal state [50] (see Eq. (32)), with logarithmic negativity [47,48] E N (t p ) = − 1 2 ln n 2 − + (n + + 1) 2 cosh 8r (A10) − (n + + 1) 4 sinh 2 8r + 4n 2 − (n + + 1) 2 cosh 2 4r , wheren ± =n 1 ±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 n − = 0, 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 2 /G 1 , i.e., the intensity of the two tones. This large mechanical entanglement is achieved when the condition e iφ 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 G for a given ∆ 0, i.e., if This dynamical scheme for the generation of continuous variable mechanical entanglement is similar to the Bogoliubov scheme proposed in Ref. [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 G and ∆ in typical situations, thereby strongly affecting the ideal Hamiltonian dynamics described here.
Appendix B: Linearization of the optomechanical dynamics with two-frequency drives The system dynamics is described by the following QLĖ where, here, differently from the description used in Sec. II, we are representing the cavity field in a reference frame rotating at the frequency ω 0 − (ω 2 − ω 1 )/2, and we have introduced the frequencies ω ± = ω 2 ± ω 1 2 .
The other parameters and operators are defined in the main text.