Stationary Gaussian Entanglement between Levitated Nanoparticles

Coherent scattering of photons is a novel mechanism of optomechanical coupling for optically levitated nanoparticles promising strong, versatile interactions with light and between nanoparticles. We show that it allows efficient deterministic generation of Gaussian entanglement between two particles in separate tweezers. A combination of red- and blue-detuned tweezers brings a mechanical Bogoliubov mode to its ground state. An additional, dispersively coupled cavity mode can reduce noise in the orthogonal mode, resulting in strong entanglement as quantified by the logarithmic negativity and verifiable with the Duan criterion for realistic experimental parameters. Such an important resource for quantum sensing and quantum simulations is pivotal for current experiments and presents an important step towards optomechanics with multiple particles in the quantum regime.

Introduction. Nonclassical states of macroscopic mechanical resonators offer a unique opportunity to investigate the boundary between classical and quantum worlds [1][2][3][4][5]. Particularly mechanical modes formed by center-of-mass motion of optically levitated particles [6][7][8] attracted attention in this regard as they do not suffer from clamping losses, allowing for efficient isolation from thermal noise. Moreover, optical trapping allows versatile modification of the potential in space and time [9][10][11]. Over the years, the field of levitated optomechanics reached remarkable level of control over the degrees of freedom down to the quantum regime [12,13], using optomechanical interactions not only for investigating the quantum-to-classical transition but also for thermodynamics [14][15][16] or force sensing [17][18][19][20]. Recently, coherent scattering of tweezer photons into an empty cavity mode arose as an interesting coupling mechanism [21], allowing efficient cooling [22,23] (even down to the quantum ground state [12]) and strong optomechanical coupling [24]. Theoretical proposals suggest to employ it for three-dimensional displacement detection [25] or for preparation of nonclassical mechanical states [26,27].
As the level of control of levitated particles progresses, systems involving multiple mechanical and electromagnetic modes become attractive for quantum sensing, simulations, and thermodynamics. These advances follow the development of cavity optomechanics with clamped mechanical resonators in which theoretical proposals and experimental demonstrations addressed generation of nonclassical correlations between fields [28][29][30][31] or between photons and phonons [32][33][34][35][36]. Entanglement between mechanical modes can also be prepared, for example, by measurements [37,38], via parametric or bluesideband driving [39][40][41][42], or using reservoir engineering [43][44][45][46]. Particularly the last strategy is attractive as it allows Gaussian entanglement to be generated deterministically in the steady state. Proposals for entanglement generation tailored for levitated systems exist as well but they either rely on weak nonlinear interactions [47] or employ photon counting [27] and are highly susceptible to thermal noise.
Here, we propose and analyze a scheme to create Gaussian two-mode squeezed entangled states of two nanoparticles deterministically in the steady state. The particles are levitated by two separate tweezers and coupled to the same optical cavity mode via coherent scattering. A suitable combination of tweezer detunings (the first tweezer is detuned to the red sideband while the second tweezer is detuned to the blue sideband; see Fig. 1) allows us to cool a collective mechanical Bogoliubov mode to its quantum ground state [43][44][45]. To combat the strong thermal noise present in the orthogonal Bogoliubov mode and open a road to entanglement generation in realistic conditions, we employ a second cavity mode coupled to both particles via the more common dispersive coupling and driven on the red sideband by an external laser [48]. Unlike existing proposals for dissipative  preparation of two-mode squeezing (which require up to four driving fields [44,46]), our scheme requires only one external drive and can generate nonclassical correlations between two mechanical modes starting from room temperature in state-of-the-art levitated systems. Our work thus presents an important step towards exploring quantum effects in the collective motion of multiple levitated particles and to the applications mentioned above.
Cooling a Bogoliubov mode via coherent scattering. We consider two particles trapped in separate tweezers and placed inside a cavity as shown in Fig. 1(a). For each particle, the tweezer defines a harmonic potential with frequency Ω j , j = 1, 2, and scattering of the tweezer photons into the cavity provides the optomechanical interaction. For the first particle (mechanical annihilation operator b 1 ), the tweezer frequency ν 1 is lower than the frequency ω 1 of cavity mode a 1 by the mechanical frequency, ω 1 = ν 1 +Ω 1 , such that scattering into the cavity mode annihilates a phonon. The other tweezer is higher than the cavity frequency by one mechanical frequency, ω 1 = ν 2 − Ω 2 , and scattering creates pairs of photons in the cavity mode and phonons in the motion of the second particle (operator b 2 ). The total Hamiltonian, under the rotating wave approximation, is given by [49] where λ j are the coupling strengths. To see how this Hamiltonian leads to steady-state two-mode squeezing, we introduce the Bogoliubov mode β 1 = ( The Hamiltonian (1) then becomes H cs = −λ eff (a † 1 β 1 + β † 1 a 1 ); this interaction cools the Bogoliubov mode β 1 to its ground state, producing two-mode squeezing between the original modes b j ; see also Fig. 1(b).
The dynamics of this system (in terms of the particle modes b j ) are described by the Langevin equationṡ where κ 1 is the decay rate of the cavity mode a 1 and γ j are the damping rates for the two mechanical modes. The input noise a 1,in has the correlations a 1,in (t)a † 1,in (t ) = δ(t − t ); the thermal noise of the mechanical modes has the correlation function b j,in (t)b † k,in (t ) = (2n j + 1)δ jk δ(t − t ), where n j k B T / Ω j is the mean thermal occupation for mechanical frequency Ω j at temperature T . To evaluate the generated entanglement, we solve for the steady-state covariance matrix associated with these dynamical equations as described in the Supplemental Material [49]. From the mechanical covariance matrix V , we then calculate the logarithmic negativity E n which is an entanglement measure for Gaussian states [50] and the state purity P = 1/ √ det V . For a direct experimental validation of the generated entanglement, we also consider the violation of the Duan criterion [51] where ∆(O) = O 2 − O 2 is the variance of the operator O and we defined the quadrature operators The violation of the Duan criterion (which is a sufficient condition for entanglement) is easier to verify experimentally as it involves only a pair of commuting operators to be measured whereas determination of the logarithmic negativity requires full tomography to obtain the whole covariance matrix V .
The results of these simulations are shown in Fig. 2(a)-(c). For these plots, we consider experimental parameters similar to the recent coherent scattering experiments [12,22]. Entanglement between the two particles can be generated for thermal noise levels surpassing n j = 10 6 as verified by the logarithmic negativity [panel (a)]; for mechanical frequency of 305 kHz, this noise level corresponds to a temperature of about 15 K. With high levels of thermal noise, the logarithmic negativity drops to zero for λ 2 → λ 1 as the coupling rate of the Bologiubov mode λ eff → 0. The EPR variance, on the other hand, drops below the classical limit only for very small amounts of thermal noise, n j 1 [panel (b)]; for n j = 10 6 (not shown), the minimum attainable EPR variance is ∼ 5 × 10 5 . This discrepancy is caused by the inherent asymmetry of the mechanical state where the first mechanical mode b 1 is actively cooled by the optomechanical interaction whereas the second mechanical mode b 2 is heated up. The Duan criterion should, in its full generality, still be able to detect entanglement when this asymmetry is taken into account. Such a verification scheme would, however, be more difficult experimentally than a direct measurement of the collective operators x 1 + x 2 and p 1 − p 2 and would require full tomography of the covariance matrix. Finally, the state purity quickly drops with increasing thermal noise as shown in panel (c). These results are understandable since we cool only the Bogoliubov mode β 1 while the orthogonal Bogoliubov mode λ eff is decoupled from the cavity field and remains in a thermal state.
Noise suppression via dispersive coupling. To make the scheme more noise-tolerant, we use optomechanical coupling to an additional cavity mode a 2 which will bring the noise in the mechanical Bogoliubov mode β 2 down. As this optical mode will be far detuned from the tweezers, we cannot use coherent scattering to couple the particles to this mode but instead employ the more conventional dispersive optomechanical coupling. We drive the mode a 2 on the red mechanical sideband (here, we assume equal mechanical frequencies, Ω 1 = Ω 2 ) (a)-(c) Generation of two-particle entanglement using coherent scattering with the Hamiltonian Hcs. We plot the logarithmic negativity (a), EPR variance (b), and mechanical-state purity (c) versus the ratio of the coupling strengths λ2/λ1 for the thermal noise levels nj = 0 (solid lines), nj = 1 (dashed), nj = 10 6 (dotted), and nj = 2 × 10 7 (dot-dashed) with the mechanical quality factor Qj = Ωj/γj = 10 8 . (d)-(f) Entanglement generation using a combination of coherent scattering and dispersive optomechanics as described by the Hamiltonian H = Hcs + H disp ; here, we consider thermal noise nj = 2 × 10 7 and mechanical quality factors Qj = 10 8 (dotted lines), Qj = 10 9 (dashed), and Qj = 10 10 (solid). The remaining parameters, equal for both sets of plots, are λ1/Ω1 = 0.35 and κ1/Ω1 = 0.4. For panels (d)-(f), the dispersive coupling strengths are g1/Ω1 = 0.01 and g2/Ω1 = 0.065 and we assume equal mechanical frequencies, Ω1 = Ω2. In panels (b), (e), the horizontal black line gives the classical limit; EPR variance smaller than this bound, ∆EPR < 1, implies that the two particles are entangled.
such that, in the rotating frame and under the rotating wave approximation, the total Hamiltonian becomes H = H cs + H disp , where H cs is given in Eq. (1) and We can express the dispersive part of the Hamiltonian as which shows that the cavity mode a 2 generates single-mode squeezing for both Bogoliubov modes [52]. Ideally, the dispersive Hamiltonian H disp should only cool the Bogoliubov mode without any additional interactions; single-mode squeezing introduces asymmetry in the steady state, making the generated entanglement undetectable by the Duan criterion. For λ 1 > λ 2 (which is necessary for dynamical stability and well-defined Bogoliubov modes) and g 1 < g 2 , the dominant part of the dispersive Hamiltonian is indeed g 2 λ 1 (a † 2 β 2 + β † 2 a 2 )/λ eff and the second cavity mode cools the Bogoliubov mode β 2 , reducing the amount of thermal noise in the mechanical steady state. The small additional squeezing contribution and a weak residual coupling to the first Bogoliubov mode β 1 slightly alleviate the steady-state noise [see also Fig. 1 We plot the resulting entanglement and mechanical state purity in Fig. 2(d)-(f). Efficient generation of entanglement-confirmed by both logarithmic negativity and EPR variance-is possible with realistic values for the thermal noise (n j 2×10 7 corresponding to thermal noise for 305 kHz mechanical modes at room temperature, T = 300 K). Particularly the drop in the EPR variance shows that the addition of the dispersive coupling helps symmetrize the mechanical steady state. This observation is further supported by the state purity which remains high for a broad range of coupling-strength ratios λ 2 /λ 1 . Finally, the considered mechanical quality factors ranging between 10 8 and 10 9 , which allow the EPR variance to drop below the classical bound of unity, have been realized in recent experiments [12,22] or can be reached with moderate improvements to the vacuum; the strength of the dispersive coupling considered here is feasible in levitated systems as well [53]. Finally, although we only show results obtained under the rotating wave approximation here, we analyse dynamics beyond this approximation in the Supplemental Material [49] and show that the finite sideband ratio still allows comparably large entanglement to be generated.
To further confirm our intuition regarding the collective Bogoliubov modes, we investigate the variation of entanglement with the dispersive coupling rates in Fig. 3. For a fixed value of the dispersive coupling g 1 , there exists an optimum coupling g 2 that simulta-neously maximizes the logarithmic negativity and purity while minimizing the EPR variance. When g 2 is too small, the cooling of the second Bogoliubov mode, ∝ g 2 λ 1 β 2 a † 2 + H.c., is weak; large g 2 , on the other hand, enhances the single-mode squeezing of the first Bogoliubov mode, ∝ g 2 λ 2 β 1 a 2 + H.c., eventually leading to instability. Moreover, as the value of g 1 starts approaching the coherent-scattering coupling (the dotted lines in Fig. 3), the generated entanglement starts to drop since the dispersive coupling begins to damp the quantum correlations generated by the coherent-scattering interaction. These results thus confirm our understanding and show that there is a broad range of dispersive coupling strengths for which entanglement between the particles can be generated.
Unequal mechanical frequencies. So far, we assumed both mechanical frequencies to be the same, Ω 1 = Ω 2 , which might be difficult to achieve in practice with nanoparticles of variable sizes. When the frequencies are different, the interactions in the coherent-scattering Hamiltonian H cs can still remain resonant when suitable tweezer frequencies are used as discussed above; the Bogoliubov mode β 1 can thus still be efficiently cooled by coherent scattering. The dispersive coupling with a single driving frequency (set on the red sideband of the second mechanical mode), on the other hand, becomes generally time-dependent and cannot efficiently cool the Bogoliubov mode β 2 . The most natural choice is then to make the interaction between a 2 and b 2 resonant as this mode is-unlike b 1 -not cooled down by the coherentscattering interaction. We thus obtain the dispersive Hamiltonian where δ 12 = Ω 1 −Ω 2 is the detuning between the two mechanical modes. To solve the resulting time-dependent Lyapunov equation in the steady state, we express the periodic dynamics in the Floquet space and solve the resulting time-independent version [49,54]. The results of this simulation (see Fig. 4) show that the detuning has negligible effect on the generated entanglement, proving that strong quantum correlations can be efficiently prepared also with particles of unequal frequencies.
If, on the other hand, the second cavity field is resonant with the first mechanical mode (with the dispersive Hamiltonian H disp = g 1 a 2 b † 1 + g 2 a 2 b † 2 e −iδ12t + H.c.), the performance of our scheme worsens, which we discuss further in the Supplemental Material [49]. The logarithmic negativity of the mechanical steady state remains high even for a large detuning between the two mechanical modes; however, as the second mechanical mode becomes off-resonant with the cavity field, the state purity decreases and the EPR variance increases. This behavior shows that this regime is not well suited for entanglement generation, further supporting our understanding that the dispersive coupling needs to efficiently remove noise from the second mechanical mode b 2 .
Discussion and outlook. For our simulations, we worked with parameters close to the recent experimental demonstrations of optomechanical cooling via coherent scattering [12,22]. The main differences are a slightly improved cavity decay rate (achievable with an improved coating of the cavity mirrors) and smaller thermal decoherence rate of the mechanics (possible by reducing the pressure in the vacuum chamber from 10 −6 mbar used in Ref. [12]). Moreover, the values of dispersive coupling considered here are also in the range available in experiments with levitated particles; coupling as strong as g/Ω = 0.1 has recently been demonstrated [53]. Finally, the effect of counterrotating terms, associated with a finite sideband ratio κ i /Ω j and coupling λ j /Ω j , does not preclude the generation of entanglement in state-ofthe-art optomechanical systems.We discuss this point in more detail in the Supplemental Material [49] where we use an extension of the Lyapunov equation into the Floquet space associated with the periodic nature of these terms [54].
In summary, we presented a scheme that generates entanglement between two levitated nanoparticles deterministically in the steady state. Coherent scattering of photons from tweezers into an empty cavity mode can be used for two-mode squeezing of the particle motion when suitable detunings between the tweezers and the cavity are used. We showed how coherent scattering from the tweezers into an empty cavity mode, together with a second, dispersively coupled cavity mode, can be used to cool two mechanical Bogoliubov modes, creating a two-mode squeezed stationary state. Crucially, our work demonstrates that it is not necessary to engineer ex-act cooling dynamics for both Bogoliubov modes which would require complicated multitone driving schemes; instead, it is sufficient to create such dynamics approximately using coherent scattering and dispersive optomechanical interaction with a single driving tone.
With levitated optomechanics entering the quantum regime, our proposal presents a viable approach towards creating complex nonclassical states of massive objects at room temperature. Our results show an attractive way towards engineering quantum dynamics and states of collective motion of multiple particles with minimal resources with possible extensions to the preparation of similar states across macroscopic distances or with more than two particles. In classical levitated optomechanics, multiparticle effects have already been observed [55]; extending the dynamics of multiple particles to the quantum regime presents an important steps forward for both basic and applied science, allowing new tests of decoherence models and schemes to measure weak forces or fields.  The total Hamiltonian of the system as depicted in Fig. 1(a) of the main text consists of the kinetic energy of the particles, the energy of the cavity fields, and the interaction between the particles and fields. As described by Gonzalez-Ballestero et al. [21], the interaction includes the trapping potential of the tweezers [given by − 1 2 αE 2 j (x j ), where E j (x j ) is the electric field of the jth tweezer at the position of the particle x j ], the usual dispersive optomechanical interaction (given by a similar expression for the electric field of the cavity mode) of the form Ga † a(b + b † ), and the coupling mediated by coherent scattering (given by the product of the tweezer and cavity fields) given by λ(a + a † )(b + b † )(e iνt + e −iνt ), where ν is the frequency of the tweezer. Moving to a frame in which the cavity mode a 1 oscillates at the frequency of the first tweezer, ν 1 , we have where ∆ 1 = ω 1 − ν 1 is the detuning between the first tweezer and the cavity resonance ω 1 and δ = ν 2 − ν 1 is the detuning between the two tweezers. Next, we set the detunings as ∆ 1 = Ω 1 , δ = Ω 1 + Ω 2 and move to the rotating frame with respect to the free oscillations of all three modes, under the rotating wave approximation, we obtain the coherent-scattering Hamiltonian Including the second cavity mode, we can rewrite Eq. (S1) as here, G j denotes the single-photon coupling rate of the dispersive optomechanical coupling and is the amplitude of an external drive at frequency ω L . Using standard linearization techniques [48] and setting the detuning ∆ 2 = ω 2 − ω L = Ω 2 , we obtain the interaction-picture Hamiltonian where g j = G j α, α is the classical steady-state amplitude of the mode a 2 , and δ 12 = Ω 1 − Ω 2 is the detuning between the two mechanical modes. Assuming equal mechanical frequencies, δ 12 = 0, the Hamiltonian (S4) simplifies to

LYAPUNOV EQUATION
The dynamics of the system are described by the Langevin equationṡ where κ j is the cavity decay rate and γ j is the mechanical amplitude damping. We also introduced the input noise operators a j,in , b j,in with the usual correlation functions with n j = k B T / Ω j denoting the average thermal occupation of the mechanical bath at temperature T . Introducing the cavity quadrature operators we can write the Langevin equations in the compact matrix forṁ where the input noise quadrature operators are defined in full analogy to the quadrature operators of the cavity and mechanical modes. The steady state of the system can now be found by solving the Lyapunov equation where V with elements V jk = r j r k + r k r j − 2 r j r j is the covariance matrix of the system's Wigner function and N = ξ(t)ξ T (t) = diag[κ 1 , κ 1 , κ 2 , κ 2 , γ 1 (2n 1 + 1), γ 1 (2n 1 + 1), γ 2 (2n 2 + 1), γ 2 (2n 2 + 1)]. (S10) The properties of the mechanical state (the logarithmic negativity, state purity, and EPR variance), are encoded in the lower right 4 × 4 block of the covariance matrix V .

DYNAMICS WITH UNEQUAL MECHANICAL FREQUENCIES
The Lyapunov equation cannot be straightforwardly used to describe a system where the two mechanical modes have different frequencies as the problem becomes time-dependent. Starting from the rotating-frame Hamiltonian (S4), we obtain the Langevin equationṡ where the vector of input noise operators is the same as before and is given by Eq. (S8c). The rotating wave approximation cannot be applied to the time-dependent part of the evolution since δ 12 = Ω 1 − Ω 2 is not necessarily larger than the cavity decay and optomechanical coupling rates.  Figure S1. Entanglement generation with particles of different mechanical frequencies. We assume that the second cavity mode a2 is driven on the red sideband of the first mechanical mode b1 as described by the Hamiltonian (S15) or drift matrix (S16). The logarithmic negativity (a), EPR variance (b), and purity (c) are plotted versus the ratio of the coupling through coherent scattering λ2/λ2 for δ12 = 0 (solid line), δ12/Ω1 = 0.2 (dashed), δ12/Ω1 = 0.4 (dot-dashed), and δ12/Ω1 = 0.8 (dotted). The remaining parameters are the same as in Fig. 4 of the main text.
To take these oscillating terms into account, we combine the Lyapunov equation with Floquet techniques to obtain an effective time-independent Lyapunov equation that includes the effect of these time-dependent terms [54]. We start by expressing the drift matrix A(t) in terms of its Fourier components defined via the expression where each of the matrices A (0) , A (n) c,s is time independent. We can now formulate the Lyapunov equation in the Floquet-space, where we introduced the Floquet-space drift matrix −δ 12 I and the drift matrix becomes (S16) The Floquet-Lyapunov approach can be applied in full analogy with the previous case; the results are shown in Fig. S1. The logarithmic negativity remains high even for detunings approaching the mechanical frequency (for which the rotating wave approximation can be applied to neglect the coupling of b 2 to the cavity mode a 2 ). As the second cavity mode is not efficiently cooled by the interaction, however, the state purity drops with increasing detuning; at the same time, the EPR variance increases, making direct verification impossible for δ 12 0.4Ω 1 .