# Generation and detection of large and robust entanglement

between two different mechanical resonators in cavity optomechanics

###### Abstract

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.

###### pacs:

42.50.Lc, 42.50.Ex, 42.50.Wk, 85.85.+j## 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 Nielsen .

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. juuls01 , while entanglement between two Josephson-junction qubits has been detected in Refs. Berkley ; Steffen . More recently, macroscopic entanglement has been demonstrated in electro-mechanical systems Palomaki : continuous variable (CV) entanglement, similar to that considered by Einstein Podolski and Rosen (EPR) epr , has been generated and detected between the position and momentum of a vibrational mode of a m-diameter Al membrane, and the quadratures of a microwave cavity field, following the theory proposal of Ref. hofer .

Entanglement between two mechanical resonators (MRs) has been instead demonstrated only at the microscopic level, in the case of two trapped ions Jost , and between two single-phonon excitations in nano-diamonds walmsley . 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 Weedbrook ; 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 marsh ; romero-isart ; ulbricht .

Many different schemes have been proposed in the literature for entangling two MRs, especially exploiting optomechanical and electromechanical devices amo ; rmp , in which the two MRs simultaneously interact with one or more electromagnetic cavity fields. Refs. PRL02 ; Peng03 ; epl considered the steady state of different systems of driven cavities: Ref. PRL02 focused on two mirrors of a ring cavity, while Ref. Peng03 assumed to drive two independent linear cavities with two-mode squeezed light transferring its entanglement to the cavity end-mirrors. Ref. epl 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 jopa , between two mechanical modes of a single movable mirror genesnjp , or in the case of two semi-transparent membranes interacting with two driven cavity modes hartmann .

A different approach for generating entangled MRs exploits conditional measurements on light modes entangled or correlated with mechanical degrees of freedom entswap ; bjorke ; mehdi1 ; woolley ; mehdi2 ; Savona . 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 JieNJP . More recent proposal applied reservoir engineering ideas poyatos ; davidovich ; zoller ; cirac ; pielawa 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 Clerk ; tan1 ; Tan ; schmidt ; WoolleyClerk ; Abdi2 ; Buchmann , eventually profiting from mechanical nonlinearities and/or parametric driving bowen-clerk ; nori .

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. Clerk ; Tan ; WoolleyClerk ; Buchmann , 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 SM . This latter scheme has been already considered in an optomechanical scenario by Kuzyk et al. Kuzyk 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 Clerk ; Tian ; Tan , 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 SM , 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 and annihilation operator interacting via the usual optomechanical interaction with two different MRs, with frequencies and and annihilation operators and respectively. The cavity mode is bichromatically driven at the two frequencies and , with the reference frequency detuned from the cavity resonance by a quantity . If we describe the cavity field in a reference frame rotating at the frequency , then the system Hamiltonian is given by

(1) |

This means that the cavity mode is simultaneously driven on the blue sideband associated with the MR with annihilation operator , and on the red sideband associated with the MR with . The nonzero detuning makes the present scheme different from the one studied in the supplementary material of Ref. Clerk which restricts to the resonant case . Our model is instead related to the scheme proposed by Kuzyk et al. Kuzyk 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, , and mechanical degrees of freedom, , are large, and one can simplify the interaction Hamiltonian at lowest order in the field fluctuations

(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 Mari , and detailed in Appendix B, approximated dynamical equations for the fluctuation operators and can be derived, in the interaction picture with respect to the Hamiltonian , 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 Eq. (79))

(3) |

The corresponding QLE including thermal noise and dissipation at rates and for the cavity and the mechanical mode respectively, are

(4) | |||

(5) | |||

(6) |

where

(7) |

are the (generally complex) linear optomechanical couplings, and and are standard input noise operators with zero mean, whose only nonzero correlation functions are , and , where is the mean thermal phonon number of the -th MR, which we assume to stay at the same environmental temperature . 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 , that we here denote with (see Appendix B). Specifically

(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 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, . In such a case, the system is stable if and only if

(9) |

This stability condition reduces to the one derived in the supplementary material of Ref. Clerk in the case . 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

(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 and real.

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

(11) | |||||

(12) |

where

(13) |

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

(14) | |||||

with the two-mode squeezing operator. The Bogoliubov mode describes the “mechanical dark mode”, which does not appear in , i.e., is decoupled from the cavity mode and therefore is a constant of motion in the absence of damping, while is the “bright” mode interacting with the cavity mode. This is equivalent to say that the dark mode 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 and . The Bogoliubov mode description has been already employed in cavity optomechanics, associated to two optical modes in Refs. Clerk ; Tian , and to two mechanical modes in Refs. Tan ; WoolleyClerk (see Appendix A for a derivation of the normal modes of the system and a study of its Hamiltonian dynamics).

### iii.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

(15) | |||

(16) | |||

(17) |

where , , are two correlated thermal noise operators whose only nonzero correlation functions are

(18) | |||||

with the effective mean thermal phonon numbers

(19) | |||||

(20) |

and the inter-mode correlation

(21) |

If a dissipative coupling term between the two Bogoliubov modes appears, which however does not have relevant effects because it is proportional to 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 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 , whose elements are is given, at the steady state, by

(22) |

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

(23) |

and

(24) |

where

(25) | |||||

(26) | |||||

(27) |

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

(28) | |||||

with

(29) |

and where now

(30) | |||||

The entanglement between modes and , measured by means of the logarithmic negativity eisert ; plenio , can be easily expressed in terms of these matrix elements as Zippilli15

(31) |

When the collective cooperativity is sufficiently large, i.e., , then 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 (), 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 occupancies and respectively. For the two MR of interest, associated with the operator and , such a state is just a two-mode squeezed thermal state marian

(32) |

where is given in Eq. (14), and

(33) |

is the density matrix of the thermal equilibrium state of a resonator with occupancy . Such a state is entangled for sufficiently large 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 , obtained from the solution of Eqs. (4)-(6). Figure 2 refers to an experimentally achievable set of parameters, s, s, s, s, and to different values of mean thermal phonon numbers , , and of the ratio . We see that remarkable values of are achieved at low temperatures, and that stationary mechanical entanglement is quite robust with respect to temperature because one has an appreciable value of even for , . 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 (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 , and a large collective cooperativity in order to significantly cool the bright Bogoliubov mode. However the corresponding optimization of the system parameters, and especially of the two couplings and , is far from being trivial. In fact, increases when , which however implies, at a fixed value of , a decreasing value of and therefore of (see Eq. (13) and Eq. (25)); moreover increasing has also the unwanted effect of increasing 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 occupancies and . At a given value of , this is obtained by taking a sufficiently large value of the associated single-mode cooperativity, , and correspondingly optimizing the value of , i.e., of . In fact, the logarithmic negativity associated with the stationary state of Eq. (32) can be evaluated in terms of the parameter

where , and can be explicitly rewritten in terms of the cooperativity as

(35) | |||||

The dependence of versus , for given values of , and , shows a maximum and then decays to zero for large (see Fig. 3 which refers to and , ). This behavior is described by a very simple approximated expression valid in the limit , with not very large , and when (corresponding to ),

(36) |

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

(37) |

given by For values of much larger or much smaller than this value, the resonators may not be entangled. When is increased to very large values , is reduced and the cooling dynamics becomes slow as compared to the standard mechanical dissipation, which takes place at rate , so that the correlations between the MRs cannot be efficiently generated. On the other hand, at small 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 and a very good estimate of the optimal value of the two-mode squeezing parameter , i.e., of , given by Eq. (37). The corresponding value of the logarithmic negativity is

(38) |

and shows that once that the ratio is optimized, the achievable stationary entanglement between the two MRs increases with increasing .

The above analysis of the stationary entanglement of the two MRs extends the results of Ref. Clerk in various directions. First of all, our model extends to the case of nonzero detuning a model discussed in the Supplementary material of Ref. Clerk . We see that a nonzero detuning has a limited effect of the dynamic of entanglement generation, providing only an effective increase of , which however becomes negligible as soon as (see Eq. (35)). Moreover, Ref. Clerk provided an explicit expression for only for the case of negligible thermal occupancies and not too large values of , while the present discussion applies for arbitrary values of , and .

## Iv Dynamical evolution in the case of equal couplings

In the special case of equal couplings , i.e., , 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. Kuzyk , 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. Kuzyk 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 milburn and by Sørensen and Mølmer SM , and first applied to an optomechanical setup by Kuzyk et al. Kuzyk . 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 , , and . One gets

(39) |

where , 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 and are two commuting conserved observables. One gets (see also Refs. Kuzyk ; milburn ; SM )

(40) | |||

(41) | |||

(42) | |||

Relevant interaction times are those when the MR dynamics decouple from that of the optical cavity, and this occurs at , , where

(43) | |||

(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

(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 and seriously affects the cavity-mediated interaction between the two MRs, as soon as becomes comparable to . Mechanical entanglement is large for large and we expect well distinct peaks for at interaction times , in the ideal parameter regime . In the more realistic regime in which , 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 s, s, s, , , and to three different values of the detuning, s (black dashed line), s (red full line), and s (blue full line). We see that an appreciable value of (even though smaller than the one achievable at the same and after the optimization of of the previous Section) is reached for a large interval of interaction times . 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. Clerk ). It is therefore interesting to study their effect by comparing the above predictions, both in the case of 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 . 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 . 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 and (that have been introduced in Eq. (II) and discussed in Appendix B), at the lowest relevant order in powers of , while for the blue ones they have been expanded at sixth order in powers of . Moreover, the green line results are found considering only the steady state solution for and , 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 .

In Fig. 5 we compare the time evolution of the entanglement evaluated with and without the time-dependent terms when . 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 and can be expanded at lowest order in . 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 and , 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 and does not provide an accurate description. We note that according to Eq. (3), in order to eliminate the fast rotating terms, the ratios have to be much larger than one. Although the dashed-blue curves are evaluated for a ratio 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 and . Specifically, when we consider either the lowest order terms and the transient dynamics, or the higher order terms and only the steady state of and , the corresponding results for the entanglement dynamics are very similar to the red lines.

In Fig. 6 we study the case of equal couplings . 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 time-dependent dynamics by calculating the Floquet exponents of the dynamical equations of the system. In fact, when and 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 Teschl ; 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 () (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 . 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. wien and then implemented in the electromechanical entanglement experiment of Ref. Palomaki , 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 , , the probe mode adiabatically follows the MR dynamics, and the output of the readout cavity is given by (see Fig. 1) wien

(46) |

with 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 .

### vi.1 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 Clerk ; Tan ; WoolleyClerk ; Kuzyk 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 . 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 .

## 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).

## Appendix A Normal modes and Hamiltonian dynamics

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

(47) |

where

(48) | |||||