Stationary Gaussian entanglement between levitated nanoparticles

We consider two particles trapped in separate tweezers and placed inside a cavity as shown in figure 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. As described by Gonzalez-Ballestero et al [23], the interaction includes the trapping potential of the tweezers (given by $-\frac{1}{2}\alpha {E}_{j}^{2}\left({x}_{j}\right)$, 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 with annihilation operator a₁) of the form ${G}_{j}{a}_{1}^{{\dagger}}{a}_{1}{x}_{j}$ (which is much weaker than coherent scattering and we neglect it in the following), and the coupling mediated by coherent scattering (given by the product of the tweezer and cavity fields) given by ${\lambda }_{j}\left({a}_{1}+{a}_{1}^{{\dagger}}\right)\left({b}_{j}+{b}_{j}^{{\dagger}}\right)\left({\text{e}}^{\text{i}{\nu }_{j}t}+{\text{e}}^{-\text{i}{\nu }_{j}t}\right)$, where ν_j is the frequency of the jth tweezer.

Moving to a frame in which the cavity mode a₁ oscillates at the frequency of the first tweezer, ν₁, we have

$$\begin{equation}H={{\Delta}}_{1}{a}_{1}^{{\dagger}}{a}_{1}+{{\Omega}}_{1}{b}_{1}^{{\dagger}}{b}_{1}+{{\Omega}}_{2}{b}_{2}^{{\dagger}}{b}_{2}-{\lambda }_{1}\left({a}_{1}+{a}_{1}^{{\dagger}}\right)\left({b}_{1}+{b}_{1}^{{\dagger}}\right)-{\lambda }_{2}\left({a}_{1}\enspace {\text{e}}^{\text{i}\delta t}+{a}_{1}^{{\dagger}}\enspace {\text{e}}^{-\text{i}\delta t}\right)\left({b}_{2}+{b}_{2}^{{\dagger}}\right),\end{equation} \tag{ 1 }$$

where Δ₁ = ω₁ − ν₁ is the detuning between the first tweezer and the cavity resonance ω₁ and δ = ν₂ − ν₁ is the detuning between the two tweezers. Next, we set the detunings as Δ₁ = Ω₁, δ = Ω₁ + Ω₂ and move to the rotating frame with respect to the free oscillations of all three modes, ${{\Delta}}_{1}{a}_{1}^{{\dagger}}{a}_{1}+{{\Omega}}_{1}{b}_{1}^{{\dagger}}{b}_{1}+{{\Omega}}_{2}{b}_{2}^{{\dagger}}{b}_{2}$; under the rotating wave approximation, we obtain the coherent-scattering Hamiltonian

$$\begin{equation}{H}_{\text{cs}}=-{\lambda }_{1}\left({a}_{1}^{{\dagger}}{b}_{1}+{b}_{1}^{{\dagger}}{a}_{1}\right)-{\lambda }_{2}\left({a}_{1}{b}_{2}+{a}_{1}^{{\dagger}}{b}_{2}^{{\dagger}}\right).\end{equation} \tag{ 2 }$$

To see how this Hamiltonian leads to steady-state two-mode squeezing, we introduce the Bogoliubov mode ${\beta }_{1}=\left({\lambda }_{1}{b}_{1}+{\lambda }_{2}{b}_{2}^{{\dagger}}\right)/{\lambda }_{\text{eff}}$, where ${\lambda }_{\text{eff}}=\sqrt{{\lambda }_{1}^{2}-{\lambda }_{2}^{2}}$ [47]. The Hamiltonian (2) then becomes ${H}_{\text{cs}}=-{\lambda }_{\text{eff}}\left({a}_{1}^{{\dagger}}{\beta }_{1}+{\beta }_{1}^{{\dagger}}{a}_{1}\right)$; this interaction cools the Bogoliubov mode β₁ to its ground state, producing two-mode squeezing between the original modes b_j (see also figure 1(b)).

The dynamics of this system (in terms of the particle modes b_j ) are described by the Langevin equations

$$\begin{equation}{\dot {a}}_{1}=\mathrm{i}\left({\lambda }_{1}{b}_{1}+{\lambda }_{2}{b}_{2}^{{\dagger}}\right)-\frac{{\kappa }_{1}}{2}{a}_{1}+\sqrt{{\kappa }_{1}}{a}_{1,\text{in}},\end{equation} \tag{ 3a }$$

$$\begin{equation}{\dot {b}}_{1}=\mathrm{i}{\lambda }_{1}{a}_{1}-\frac{{\gamma }_{1}}{2}{b}_{1}+\sqrt{{\gamma }_{1}}{b}_{1,\text{in}},\end{equation} \tag{ 3b }$$

$$\begin{equation}{\dot {b}}_{2}=\mathrm{i}{\lambda }_{2}{a}_{1}^{{\dagger}}-\frac{{\gamma }_{2}}{2}{b}_{2}+\sqrt{{\gamma }_{2}}{b}_{2,\text{in}},\end{equation} \tag{ 3c }$$

where κ₁ is the decay rate of the cavity mode a₁ and γ_j are the damping rates for the two mechanical modes. The input noise a_1,in has the correlations $\langle {a}_{1,\text{in}}\left(t\right){a}_{1,\text{in}}^{{\dagger}}\left({t}^{\prime }\right)\rangle =\delta \left(t-{t}^{\prime }\right)$; the thermal noise of the mechanical modes has the correlation function $\langle {b}_{j,\text{in}}\left(t\right){b}_{k,\text{in}}^{{\dagger}}\left({t}^{\prime }\right)\rangle =\left(2{n}_{j}+1\right){\delta }_{jk}\delta \left(t-{t}^{\prime }\right)$, where n_j ≃ k_B T/ℏΩ_j is the mean thermal occupation for mechanical frequency Ω_j at temperature T. In terms of the cavity quadrature operators ${X}_{1}=\left({a}_{1}+{a}_{1}^{{\dagger}}\right)/\sqrt{2},{Y}_{1}=-\mathrm{i}\left({a}_{1}-{a}_{1}^{{\dagger}}\right)/\sqrt{2}$ and canonical mechanical operators ${x}_{j}=\left({b}_{j}+{b}_{j}^{{\dagger}}\right)/\sqrt{2},{p}_{j}=-\mathrm{i}\left({b}_{j}-{b}_{j}^{{\dagger}}\right)/\sqrt{2}$, the dynamics can be expressed in the form

$$\begin{equation}\dot {r}=Ar+\xi ,\end{equation} \tag{ 4a }$$

$$\begin{equation}r={\left({X}_{1},{Y}_{1},{x}_{1},{p}_{1},{x}_{2},{p}_{2}\right)}^{T},\end{equation} \tag{ 4b }$$

$$\begin{equation}\xi ={\left(\sqrt{{\kappa }_{1}}{X}_{1,\text{in}},\sqrt{{\kappa }_{1}}{Y}_{1,\text{in}},\sqrt{{\gamma }_{1}}{x}_{1,\text{in}},\sqrt{{\gamma }_{1}}{p}_{1,\text{in}},\sqrt{{\gamma }_{2}}{x}_{2,\text{in}},\sqrt{{\gamma }_{2}}{p}_{2,\text{in}}\right)}^{T}\end{equation} \tag{ 4c }$$

$$\begin{equation}A=\left(\begin{matrix}\hfill -\frac{1}{2}{\kappa }_{1}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -{\lambda }_{1}\hfill & \hfill 0\hfill & \hfill {\lambda }_{2}\hfill \\ \hfill 0\hfill & \hfill -\frac{1}{2}{\kappa }_{1}\hfill & \hfill {\lambda }_{1}\hfill & \hfill 0\hfill & \hfill {\lambda }_{2}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -{\lambda }_{1}\hfill & \hfill -\frac{1}{2}{\gamma }_{1}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill {\lambda }_{1}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -\frac{1}{2}{\gamma }_{1}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\lambda }_{2}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -\frac{1}{2}{\gamma }_{2}\hfill & \hfill 0\hfill \\ \hfill {\lambda }_{2}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -\frac{1}{2}{\gamma }_{2}\hfill \end{matrix}\right),\end{equation} \tag{ 4d }$$

where the input noise quadrature operators are defined in full analogy to the quadrature operators of the cavity and mechanical modes.

To evaluate the generated entanglement, we solve for the steady-state covariance matrix associated with these dynamical equations as described in the appendix. From the mechanical covariance matrix V, we then calculate the logarithmic negativity E_n which is an entanglement measure for Gaussian states [52] and the state purity $P=1/\sqrt{\mathrm{det}\enspace V}$. For a direct experimental validation of the generated entanglement, we also consider the violation of the Duan criterion [53]

$$\begin{equation}{{\Delta}}_{\text{EPR}}=\frac{1}{2}\left[{\Delta}\left({x}_{1}+{x}_{2}\right)+{\Delta}\left({p}_{1}-{p}_{2}\right)\right]{\geqslant}1,\end{equation} \tag{ 5 }$$

where Δ(O) = ⟨O²⟩ − ⟨O⟩² is the variance of the operator O. 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. To further highlight the difference between our results and previous work [29], we also calculate the noise reduction factor [54] which is defined as the variance of the phonon-number difference normalised by the total number of excitations,

$$\begin{equation}\text{NRF}=\frac{\langle {\left({n}_{1}-{n}_{2}\right)}^{2}\rangle -{\langle \left({n}_{1}-{n}_{2}\right)\rangle }^{2}}{\langle {n}_{1}+{n}_{2}\rangle },\end{equation} \tag{ 6 }$$

where the phonon number ${n}_{j}={b}_{j}^{{\dagger}}{b}_{j}$. Noise reduction factor smaller than unity is a clear signature of sub-Poissonian phonon–phonon correlations (and it vanishes completely for pure two-mode squeezing). Conversely, single-phonon entangled states proposed in [29] exhibit phonon anti-correlation between particles with NRF = 1.

The results of these simulations are shown in figures 2(a)–(d). For these plots, we consider experimental parameters similar to the recent coherent scattering experiments [12, 24]. Entanglement between the two particles can be generated for thermal noise levels surpassing n_j = 10⁶ 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 λ₂ → λ₁ 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⁶ (not shown), the minimum attainable EPR variance is ∼5 × 10⁵. This discrepancy is caused by the inherent asymmetry of the mechanical state where the first mechanical mode b₁ is actively cooled by the optomechanical interaction whereas the second mechanical mode b₂ 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₁ + x₂ and p₁ − p₂ and would require full tomography of the covariance matrix. This asymmetry is further highlighted by the noise reduction factor which remains smaller than unity only for extremely small amounts of thermal noise (panel (c)). Finally, the state purity quickly drops with increasing thermal noise as shown in panel (d). These results are understandable since we cool only the Bogoliubov mode β₁ while the orthogonal Bogoliubov mode ${\beta }_{2}=\left({\lambda }_{2}{b}_{1}^{{\dagger}}+{\lambda }_{1}{b}_{2}\right)/\sqrt{{\lambda }_{\text{eff}}}$ is decoupled from the cavity field and remains in a thermal state.

Zoom In Zoom Out Reset image size

To make the scheme more noise-tolerant, we use optomechanical coupling to an additional cavity mode a₂ which will bring the noise in the mechanical Bogoliubov mode β₂ 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₂ on the red mechanical sideband (here, we assume equal mechanical frequencies, Ω₁ = Ω₂) 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 (2) and

$$\begin{equation}{H}_{\text{disp}}={g}_{1}\left({a}_{2}^{{\dagger}}{b}_{1}+{b}_{1}^{{\dagger}}{a}_{2}\right)+{g}_{2}\left({a}_{2}^{{\dagger}}{b}_{2}+{b}_{2}^{{\dagger}}{a}_{2}\right).\end{equation} \tag{ 7 }$$

We can express the dispersive part of the Hamiltonian as ${H}_{\text{disp}}=\left({g}_{1}{\lambda }_{1}{\beta }_{1}-{g}_{2}{\lambda }_{2}{\beta }_{1}^{{\dagger}}\right){a}_{2}^{{\dagger}}/{\lambda }_{\text{eff}}+\left({g}_{2}{\lambda }_{1}{\beta }_{2}-{g}_{1}{\lambda }_{2}{\beta }_{2}^{{\dagger}}\right){a}_{2}^{{\dagger}}/{\lambda }_{\text{eff}}+\text{H.c.}$ which shows that the cavity mode a₂ generates single-mode squeezing for both Bogoliubov modes [55]. 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 λ₁ > λ₂ (which is necessary for dynamical stability and well-defined Bogoliubov modes) and g₁ < g₂, the dominant part of the dispersive Hamiltonian is indeed ${g}_{2}{\lambda }_{1}\left({a}_{2}^{{\dagger}}{\beta }_{2}+{\beta }_{2}^{{\dagger}}{a}_{2}\right)/{\lambda }_{\text{eff}}$ and the second cavity mode cools the Bogoliubov mode β₂, 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 β₁ slightly alleviate the steady-state noise (see also figure 1(b)).

We plot the resulting entanglement and mechanical state purity in figures 2(e)–(h). Efficient generation of entanglement—confirmed by logarithmic negativity, EPR variance, as well as the noise reduction factor—is possible with realistic values of thermal noise (n_j ≃ 2 × 10⁷ corresponding to thermal noise for 305 kHz mechanical modes at room temperature, T = 300 K). Particularly the drop in the EPR variance and noise reduction factor shows that the addition of the dispersive coupling helps to symmetrise the mechanical steady state. This observation is further supported by the state purity which remains high for a broad range of coupling-strength ratios λ₂/λ₁. Finally, the considered mechanical quality factors ranging between 10⁸ and 10⁹, which allow the EPR variance to drop below the classical bound of unity, have been realised in recent experiments [12, 24] 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 [56].

To further confirm our intuition regarding the collective Bogoliubov modes, we investigate the variation of entanglement with the dispersive coupling rates in figure 3. For a fixed value of the dispersive coupling g₁, there exists an optimum coupling g₂ that simultaneously maximises the logarithmic negativity and purity while minimizing the EPR variance. When g₂ is too small, the cooling of the second Bogoliubov mode, $\propto {g}_{2}{\lambda }_{1}{\beta }_{2}{a}_{2}^{{\dagger}}+\text{H.c.}$, is weak; large g₂, on the other hand, enhances the single-mode squeezing of the first Bogoliubov mode, ∝g₂ λ₂ β₁ a₂ + H.c., eventually leading to instability. Moreover, as the value of g₁ starts approaching the coherent-scattering coupling (the dotted lines in figure 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.

Zoom In Zoom Out Reset image size

So far we have shown the results obtained under the rotating wave approximation. To provide a more detailed analysis of our proposal under realistic experimental conditions, we now analyse two situations where the rotating wave approximation cannot be applied. First, we study the effect of counterrotating terms omitted from the Hamiltonian (2) and show that entanglement can still be generated when realistic values for the sideband ratio κ_j /Ω₁ and coupling λ_j /Ω₁ are included. Next, we consider a scenario where the two mechanical frequencies are different, Ω₁ ≠ Ω₂, and demonstrate that our scheme is capable of generating strong entanglement in this case as well.

To analyse dynamics beyond the rotating wave approximation, we start from the Hamiltonian (1) appended with the dispersive part of the interaction, ${H}_{\text{disp}}={{\Delta}}_{2}{a}_{2}^{{\dagger}}{a}_{2}+{g}_{1}\left({a}_{2}+{a}_{2}^{{\dagger}}\right)\left({b}_{1}+{b}_{1}^{{\dagger}}\right)+{g}_{2}\left({a}_{2}+{a}_{2}^{{\dagger}}\right)\left({b}_{2}+{b}_{2}^{{\dagger}}\right)$, and move to the rotating frame with respect to the free oscillations, ${H}_{0}={{\Delta}}_{1}{a}_{1}^{{\dagger}}{a}_{1}+{{\Delta}}_{2}{a}_{2}^{{\dagger}}{a}_{2}+{{\Omega}}_{1}{b}_{1}^{{\dagger}}{b}_{1}+{{\Omega}}_{2}{b}_{2}^{{\dagger}}{b}_{2}$. Assuming again Ω₁ = Ω₂, we obtain the full interaction-picture Hamiltonian

$$\begin{equation}\begin{aligned}\hfill H& =-{\lambda }_{1}\left({a}_{1}^{{\dagger}}{b}_{1}+{b}_{1}^{{\dagger}}{a}_{1}\right)-{\lambda }_{2}\left({a}_{1}{b}_{2}+{a}_{1}^{{\dagger}}{b}_{2}^{{\dagger}}\right)-{\lambda }_{1}\left({a}_{1}{b}_{1}\enspace {\mathrm{e}}^{-2\mathrm{i}{{\Omega}}_{1}t}+{a}_{1}^{{\dagger}}{b}_{1}^{{\dagger}}\enspace {\mathrm{e}}^{2\mathrm{i}{{\Omega}}_{1}t}\right)-{\lambda }_{2}\left({a}_{1}{b}_{2}^{{\dagger}}\enspace {\mathrm{e}}^{2\mathrm{i}{{\Omega}}_{1}t}+{a}_{1}^{{\dagger}}{b}_{2}\enspace {\mathrm{e}}^{-2\mathrm{i}{{\Omega}}_{1}t}\right)\hfill \\ \hfill & \quad +{g}_{1}\left({a}_{2}^{{\dagger}}{b}_{1}+{b}_{1}^{{\dagger}}{a}_{2}\right)+{g}_{2}\left({a}_{2}^{{\dagger}}{b}_{2}+{b}_{2}^{{\dagger}}{a}_{2}\right)+{g}_{1}\left({a}_{2}{b}_{1}\enspace {\mathrm{e}}^{-2\mathrm{i}{{\Omega}}_{1}t}+{a}_{2}^{{\dagger}}{b}_{1}^{{\dagger}}\enspace {\mathrm{e}}^{2\mathrm{i}{{\Omega}}_{1}t}\right)+{g}_{2}\left({a}_{2}{b}_{2}\enspace {\mathrm{e}}^{-2\mathrm{i}{{\Omega}}_{1}t}+{a}_{2}^{{\dagger}}{b}_{2}^{{\dagger}}\enspace {\mathrm{e}}^{2\mathrm{i}{{\Omega}}_{1}t}\right),\hfill \end{aligned}\end{equation} \tag{ 8 }$$

from which we express the time-dependent Langevin equations for the quadrature operators

$$\begin{equation}\dot {r}=A\left(t\right)r+\xi ,\end{equation} \tag{ 9a }$$

$$\begin{equation}A\left(t\right)=\left(\begin{matrix}\hfill -\frac{1}{2}{\kappa }_{1}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -{\lambda }_{1}s\left(t\right)\hfill & \hfill -{\lambda }_{1}{c}_{-}\left(t\right)\hfill & \hfill {\lambda }_{2}s\left(t\right)\hfill & \hfill {\lambda }_{2}{c}_{-}\left(t\right)\hfill \\ \hfill 0\hfill & \hfill -\frac{1}{2}{\kappa }_{1}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {\lambda }_{1}{c}_{+}\left(t\right)\hfill & \hfill {\lambda }_{1}s\left(t\right)\hfill & \hfill {\lambda }_{2}{c}_{+}\left(t\right)\hfill & \hfill {\lambda }_{2}s\left(t\right)\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -\frac{1}{2}{\kappa }_{2}\hfill & \hfill 0\hfill & \hfill {g}_{1}s\left(t\right)\hfill & \hfill {g}_{1}{c}_{-}\left(t\right)\hfill & \hfill {g}_{2}s\left(t\right)\hfill & \hfill {g}_{2}{c}_{-}\left(t\right)\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -\frac{1}{2}{\kappa }_{2}\hfill & \hfill -{g}_{1}{c}_{+}\left(t\right)\hfill & \hfill -{g}_{1}s\left(t\right)\hfill & \hfill -{g}_{2}{c}_{+}\left(t\right)\hfill & \hfill -{g}_{2}s\left(t\right)\hfill \\ \hfill -{\lambda }_{1}s\left(t\right)\hfill & \hfill -{\lambda }_{1}{c}_{-}\left(t\right)\hfill & \hfill {g}_{1}s\left(t\right)\hfill & \hfill {g}_{1}{c}_{-}\left(t\right)\hfill & \hfill -\frac{1}{2}{\gamma }_{1}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill {\lambda }_{1}{c}_{+}\left(t\right)\hfill & \hfill {\lambda }_{1}s\left(t\right)\hfill & \hfill -{g}_{1}{c}_{+}\left(t\right)\hfill & \hfill -{g}_{1}s\left(t\right)\hfill & \hfill 0\hfill & \hfill -\frac{1}{2}{\gamma }_{1}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill -{\lambda }_{2}s\left(t\right)\hfill & \hfill {\lambda }_{2}{c}_{-}\left(t\right)\hfill & \hfill {g}_{2}s\left(t\right)\hfill & \hfill {g}_{2}{c}_{-}\left(t\right)\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -\frac{1}{2}{\gamma }_{2}\hfill & \hfill 0\hfill & \hfill \hfill \\ \hfill {\lambda }_{2}{c}_{+}\left(t\right)\hfill & \hfill -{\lambda }_{2}s\left(t\right)\hfill & \hfill -{g}_{2}{c}_{+}\left(t\right)\hfill & \hfill -{g}_{2}s\left(t\right)\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -\frac{1}{2}{\gamma }_{2}\hfill \end{matrix}\right),\end{equation} \tag{ 9b }$$

where c_±(t) = 1 ± cos(2Ω₁ t), s(t) = sin(2Ω₁ t). To find the steady state in this case, we express the drift matrix A(t) in the Floquet space [57]. Using the approach outlined in the appendix, we can thus turn the time-dependent problem into a time-independent one for which the corresponding Lyapunov equation can be solved to obtain the steady state covariance matrix.

Results of the corresponding simulations are shown in figure 4. The finite sideband ratio (κ_j /Ω₁ ≃ 0.4 for the experiment in [12]) and strong coupling (λ₁/Ω₁ ≃ 0.35) allow generation of entanglement even when the effect of counterrotating terms is taken into account. Crucially, the noise due to the counterrotating terms is distributed evenly between the two mechanical modes, allowing the presence of entanglement to be verified from the symmetric EPR variance.

Zoom In Zoom Out Reset image size

Another situation where the rotating wave approximation cannot be applied arises when the two mechanical frequencies are different, Ω₁ ≠ Ω₂, which might be the case for two nanoparticles of different sizes. In this case, 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 β₁ 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 β₂.

The most natural choice is then to make the interaction between the modes a₂ and b₂ resonant as this mechanical mode is—unlike b₁—not cooled down by the coherent-scattering interaction. We thus obtain the dispersive Hamiltonian

$$\begin{equation}{H}_{\text{disp}}={g}_{1}{a}_{2}{b}_{1}^{{\dagger}}\enspace {\text{e}}^{\text{i}{\delta }_{12}t}+{g}_{2}{a}_{2}{b}_{2}^{{\dagger}}+\text{H.c.},\end{equation} \tag{ 10 }$$

where δ₁₂ = Ω₁ − Ω₂ is the detuning between the two mechanical modes and we applied the rotating wave approximation to neglect terms oscillating at 2Ω_1,2 and Ω₁ + Ω₂. To solve the resulting time-dependent Lyapunov equation in the steady state, we again express the periodic dynamics in the Floquet space and solve the resulting time-independent version as described in the appendix. The results of this simulation in figures 5(a)–(c) show that the detuning has no observable effect on the generated entanglement, proving that strong quantum correlations can be efficiently prepared also with particles of unequal frequencies.

Zoom In Zoom Out Reset image size

If, on the other hand, the second cavity field is resonant with the first mechanical mode as described by the Hamiltonian

$$\begin{equation}{H}_{\text{disp}}={g}_{1}{a}_{2}{b}_{1}^{{\dagger}}+{g}_{2}{a}_{2}{b}_{2}^{{\dagger}}\enspace {\text{e}}^{-\text{i}{\delta }_{12}t}+\text{H.c.},\end{equation} \tag{ 11 }$$

the performance of our scheme worsens which we show in figures 5(d)–(f). 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 behaviour 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₂.

For our simulations, we worked with parameters close to the recent experimental demonstrations of optomechanical cooling via coherent scattering [12, 24]. 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⁻⁶ mbar used in reference [12]). Note, however, that our model does not assume decoherence merely by gas damping. In our effective, phase-insensitive decoherence we include a broad range of effects. The crucial parameter in our simulations is the thermal decoherence rate γn_j which ought to be compared with experimentally available heating rates. Admittedly, our heating rate γ_j n_j /2π ≃ 6 kHz (corresponding to Q_j = Ω/γ_j = 10⁹) is smaller than the total heating rate reported in reference [12]; nevertheless, it corresponds to the recoil heating rate reported there.

Moreover, the values of dispersive coupling considered here are also within reach of current experiments. Coupling approaching g/Ω ≃ 0.1 has recently been demonstrated, albeit with a smaller mechanical frequency (of about 163 kHz) and coupling rate (of about 14.4 kHz) [56]. Given the broad range of experimental parameters over which nonclassical two-mode squeezed states can be prepared (see figure 3), our proposed scheme is feasible.

Importantly, decoherence processes introduced by the cavity drive do not pose a limitation in our proposal. First, photon recoil, proportional to the light intensity, is much weaker than the recoil caused by the trapping beam. The intensity of the tweezer is about 8.2 × 10⁷ W cm⁻² (for power P = 0.4 W and waist w = 0.7 μm [12]) whereas the cavity-field intensity is merely 1.9 × 10⁶ W cm⁻² (P = 30 W, w = 40 μm; this value corresponds to n_phot = 10¹⁰ photons needed to reach coupling of 30 kHz from single-photon coupling of 0.3 Hz [56]). The heating rate associated with the cavity photon recoil can therefore be expected to be of the order of 100 Hz.

A bigger problem associated with the cavity drive is the laser phase noise; the corresponding heating rate can be expressed as

$$\begin{equation*}{{\Gamma}}_{\text{phase}}=\frac{4{g}^{2}{n}_{\text{phot}}}{{\kappa }^{2}}{S}_{\dot {\phi }}\left({\Omega}\right),\end{equation*}$$

where ${S}_{\dot {\phi }}\left({\Omega}\right)$ is the frequency noise spectrum. With typical values of 0.1 Hz² Hz⁻¹, we get a heating rate of the order of 100 MHz. The laser phase noise is, however, only a technical limitation and can be further reduced with the help of a filtering cavity [58]; to obtain heating rate comparable to the photon recoil from the tweezer, a filtering cavity with a linewidth of about 40 kHz is sufficient.

The generated entanglement can be most easily detected via the EPR variance (5) which requires only estimating the variance of the sum of the mechanical positions, Δ(x₁ + x₂), and of the difference of mechanical momenta, Δ(p₁ − p₂), and not tomography of the full covariance matrix. The best strategy in terms of measurement precision would employ a third cavity mode (such that the above described interactions still actively stabilise the entangled state) using backaction evading measurements of a single joint mechanical quadrature [59]. Such a measurement allows, in principle, sensitivity below the standard quantum limit as measurement backaction affects the two orthogonal collective quadratures x₁ − x₂, p₁ + p₂; since the two measured quadratures commute, they can be measured simultaneously (or subsequently) with shot-noise limited precision.

The measurement requires both particles to be coupled to the readout cavity mode with the same rate. This arrangement ensures that the collective mechanical quadrature coupled to the cavity field contains both particles with the same weight. As it might be difficult to achieve such coupling via radiation pressure interaction (which would introduce additional photon recoil and laser phase noise), the readout could be performed via coherent scattering as well. Using a pair of weak additional tweezer beams for each particle—one on the blue sideband, one on the red—would lead to scattering of both Stokes and anti-Stokes photons into the cavity, analogous to two-tone driving in conventional optomechanics [59]. Such a setting would also have the advantage that, for each particle, the mechanical quadrature coupled to the cavity field would be set by the relative phase between the two tweezers; the measured quadratures for the two particles can thus be set independently. This readout therefore allows also measurement of the orthogonal collective quadratures x₁ − x₂, p₁ + p₂, as well as single-particle quadratures, and thus enables full quantum tomography of the mechanical state which can be used to estimate the noise reduction factor (6) or the logarithmic negativity.