Introduction

Entanglement is a crucial resource for quantum information processing (QIP)1. As such, the ability to entangle fields of arbitrary wavelength will be important for linking nodes in heterogeneous QIP networks2. Mechanical oscillators are uniquely poised in their ability to create such links, thanks to the frequency-independence of the radiation pressure interaction3,4 and their ability to couple to various solid-state qubits5,6,7,8,9,10. The ability to entangle two radiation fields via a common mechanical interaction was outlined 20 years ago11,12,13, and the intervening decades have seen the development of optomechanical devices14 which are robustly quantum mechanical and routinely integrated into hybrid systems.

Recently, mechanically mediated entanglement has been reported between propagating microwave fields15 as well as two superconducting qubits16. In both cases, the entanglement remained confined to the dilution refrigerator in which it was created. Here, we create mechanically entangled optical fields for the first time, spanning up to 100 nm in wavelength, by utilizing an extremely coherent mechanical platform (Fig. 1). Moreoever, while the entangling mechanical interface resides at cryogenic temperatures, it is compatible with highly efficient light extraction and collection, such that we can directly measure the entanglement at room temperature, without noise subtraction or other indirect inference. This in turn means that the entangled optical fields could easily be distributed for QIP applications.

Results

Theoretical model

We consider two propagating optical fields (labeled by j = AB), from which one can identify a pair of temporal modes with quadratures \({\hat{X}}_{j},{\hat{Y}}_{j}\). We take the variance of these modes to be 1/2 in their ground state. From these modes, one can construct joint Einstein–Podolsky–Rosen (EPR)-type variables \({\hat{X}_{\pm }}={\hat{X}_{A}}\pm {\hat{X}_{B}}\) and \({\hat{Y}_{\pm }}={\hat{Y}_{A}}\pm {\hat{Y}_{B}}\), which form the basis for various entanglement criteria17,18. We adopt the common Duan–Giedke–Cirac–Zoller (DGCZ) criterion17 for the inseparability \({\mathcal{I}}\), which states that the two modes are inseparable if their variances (V) satisfy

$${\mathcal{I}}\equiv \frac{V({\hat{X}}_{+})+V({\hat{Y}}_{-})}{2} \, < \, 1.$$
(1)

To further quantify this entanglement, one can utilize the system’s covariance matrix, σ, which fully characterizes the correlations between various quadratures. From this matrix, it is straightforward to calculate the symplectic eigenvalues of its partial transpose, \({\tilde{\nu }_{\pm }}\!\)19. These eigenvalues offer a condition for separability (\(2{\tilde{\nu }}_{-}\ge 1\)), as well as a tool to calculate a common measure of entanglement, the logarithmic negativity \({E}_{N}=\max \left[0,-{\mathrm{log}\,}_{2}\ 2{\tilde{\nu }}_{-}\right]\). We note that 2\({\tilde{\nu }}_{-}\) also corresponds to the minimum value of \({\mathcal{I}}\) possible when optimizing over local operations on either subsystem (e.g., squeezing and rotation)20. Thus, 2\({\tilde{\nu }}_{-}\) serves as a lower bound for any DGCZ measurement.

Fig. 1: Experimental setup.
figure 1

Two lasers (red and blue) simultaneously drive an optomechanical cavity, kept in a helium flow cryostat. The inset shows the structure of the soft-clamped mechanical resonator (Si3N4 in white, holes in black). Exiting the cavity, the optical fields possess nonlocal correlations, illustrated by the squeezed phase space ellipses. After the cavity, the two lasers are physically separated and detected simultaneously by balanced homodyne detectors, with local oscillators locked at phases θAθB. The top of the figure shows a frequency diagram of the relevant optical modes. The two cavity drives are shown in black, with scattered mechanical sidebands of laser A and B shown in blue and red, respectively. The sideband quadrature modes considered in the paper correspond to combinations of both scattered sidebands, as indicated by the blue and red shaded areas.

In an optomechanical setting, in which a mechanical resonator (unitless position \(\hat{q}\), momentum \(\hat{p}\)) is linearly coupled to two independent, resonantly driven cavity modes (amplitude \({\hat{X}}_{j}^{{\rm{cav}}}\), phase \({\hat{Y}}_{j}^{{\rm{cav}}}\), j = AB), the interaction Hamiltonian can be written as: \({\hat{H}}_{{\rm{int}}}=-{\sum }_{j}2\hslash {g}_{j}{\hat{X}}_{j}^{{\rm{cav}}}\hat{q}\), where is the reduced Planck constant, and gj are the field-enhanced optomechanical coupling rates. We consider an unresolved-sideband system, in which the cavity decay rates, κj are much faster than the mechanical frequency Ωm and mechanical energy dissipation rate Γm. We also assume that the two cavity modes are driven symmetrically by coherent states, such that their induced quantum backaction rates, \({\Gamma }_{j}^{{\rm{qba}}}=4{g}_{j}^{2}/{\kappa }_{j}\) are equal: \({\Gamma }_{A}^{{\rm{qba}}}={\Gamma }_{B}^{{\rm{qba}}}\equiv {\Gamma }_{{\rm{qba}}}\). The quadratic interaction preserves the Gaussianity of the state. The following equations of motion link the input and output optical fields

$${\hat{X}}_{j}^{{\rm{out}}}(t)=-{\hat{X}}_{j}^{{\rm{in}}}(t),$$
(2a)
$${\hat{Y}}_{j}^{{\rm{out}}}(t)= -{\hat{Y}}_{j}^{{\rm{in}}}(t)-2\sqrt{{\Gamma }_{{\rm{qba}}}}{\chi }_{{m}}(t)* [\sqrt{2{\Gamma }_{m}}{\hat{P}}_{{\rm{in}}}(t) \\ +{\sum _{i = A,B}}2\sqrt{{\Gamma }_{{\rm{qba}}}}{\hat{X}}_{i}^{{\rm{in}}}(t)],$$
(2b)

where \({\hat{X}}_{j}^{{\rm{in}}}\) and \({\hat{Y}}_{j}^{{\rm{in}}}\) are input vacuum noise quadratures, χm is the mechanical susceptibility, Γm is the mechanical energy dissipation rate, * indicates convolution, and \({\hat{P}}_{{\rm{in}}}\) is the mechanical thermal noise operator.

In Eq. (2), we see that the quantum amplitude fluctuations of each laser drive the mechanical system, whose motion is then imprinted on the optical phase. This is the same mechanism that drives ponderomotive squeezing of a single laser21, but in this case there are also cross-correlations between the lasers. More insight can be had by moving to the joint mode basis (see Supplementary Note 2), where one finds that the system decouples into a sum mode undergoing ponderomotive squeezing, and a difference mode which remains dark to all mechanical dynamics. It is this squeezing of a joint (nonlocal) mode which results in the “ponderomotive entanglement” we study here. We note that Eq. (2) also closely mirror those describing four-mode squeezing based on the Kerr nonlinearity in glass22, in which the response function χm is effectively instantaneous. Both approaches are members of a broader class of settings that enable, in principle, quantum-non-demolition measurements of light23,24.

Homodyne detection allows measurement of optical quadratures in a rotated basis defined by the local oscillator phase, θj. Filtering the homodyne signal at frequency Ω with a mode function h(t) yields the sideband quadratures of a particular temporal mode at time t:

$${\hat{X}}_{j}^{{\theta }_{j}}(t)={\int _{-\infty }^{t}}ds\ \cos (\Omega s)h(t-s)\left({\hat{X}}_{j}^{{\rm{out}}}(s)\cos ({\theta }_{j})+{\hat{Y}}_{j}^{{\rm{out}}}(s)\sin ({\theta }_{j})\right),$$
(3a)
$${\hat{Y}}_{j}^{{\theta }_{j}}(t)={\hat{X}}_{j}^{{\theta }_{j}+\pi /2}(t),$$
(3b)

where θj is the homodyne angle. Note that the quadratures \({\hat{X}}_{j}^{{\theta }_{j}}\) available in a homodyne detector contain a pair of sidebands, symmetric to the carrier, as illustrated at the top of Fig. 1. As there are in total four optical modes involved, correlations between such modes are sometimes called four-mode-squeezing (see Supplementary Note 3), in contrast to the entangled microwave modes recently analyzed in a heterodyne scheme15. In the following model, we consider the limit of long filter times, in which h effectively selects a single Fourier component20. Furthermore, since the system is stationary, we drop the time argument t and focus on the ensemble statistics of these modes.

Within this model, one can obtain a simple expression for the DGCZ inseparability criterion (see Supplementary Note 1 for detail)

$${{\mathcal{I}}}_{\Theta ,\Omega }^{{\rm{ideal}}}= \,\,1+8{\Gamma }_{{\rm{qba}}}| {\chi }_{{m}}(\Omega ){| }^{2}{\Gamma }_{{\rm{dec}}}\left(1+\cos (2\Theta )\right)\\ -4{\Gamma }_{{\rm{qba}}}{\rm{Re}}\left[{\chi }_{{m}}(\Omega )\right]\sin (2\Theta ),$$
(4)

where Θ ≡ (θA + θB)/2. The first term is the contribution from shot noise at the detectors. The second term is the contribution from mechanical motion, where the total decoherence rate \({\Gamma }_{{\rm{dec}}}=2{\Gamma }_{{\rm{qba}}}+{\Gamma }_{m}({\bar{n}}_{{\rm{th}}}+1/2)\) includes both quantum backaction sources and thermal motion. The third term corresponds to correlation between two beams, again in close analogy to ponderomotive squeezing21. In practice, there is always optical loss, which admits vacuum noise that degrades the detected correlations. This is described by a collection efficiency ηA = ηB ≡ η < 1, with which the inseparability of the detected optical states becomes \({{\mathcal{I}}}_{\Theta ,\Omega }=\eta {{\mathcal{I}}}_{\Theta ,\Omega }^{{\rm{ideal}}}+(1-\eta )\). By defining a combined measurement efficiency ηmeas = 2ηΓqba/Γdec, one can show (see Supplementary Note 1) that the minimum of \({{\mathcal{I}}}_{\Theta ,\Omega }\) is given by 1 − ηmeas/2. By further calculating the full covariance matrix for this toy model (see Supplementary Note 1), one can show that \(\min \{2{\tilde{\nu }}_{-}\}=\sqrt{1-{\eta }_{{\rm{meas}}}}\), that is, the system can generate arbitrarily strong entanglement as ηmeas → 1.

Experimental setup

In practice, the optical fields become entangled via their shared interaction with a 3.6 mm × 3.6 mm × 20 nm soft-clamped Si3N4 membrane25. The vibrational mode of central defect has a frequency of Ωm = 2π × 1.139 MHz, and a quality factor Q = 1.04 × 109 at a temperature of 10 K, which corresponds to a mechanical linewidth of Γm = 2π × 1.10 mHz.

As illustrated in Fig. 1, the membrane is inserted in the middle of an optical cavity26,27, addressed by two lasers with wavelength  ~796 nm. These lasers are orthogonally polarized and populate the cavity in two different longitudinal modes separated by  ~0.3 THz, with linewidths of κA = 2π × 13.3 MHz and κB = 2π × 12.6 MHz. With this setup, we achieve \({\Gamma }_{A}^{{\rm{qba}}}\approx 2\pi \times 1.35\) kHz and \({\Gamma }_{B}^{{\rm{qba}}}\approx 2\pi \times 0.89\) kHz, which easily exceed the thermal decoherence rate \({\Gamma }_{{\rm{m}}}{\bar{n}}_{{\rm{th}}}\approx 2\pi \times 0.20\) kHz. We measure the optical quadratures of the cavity output fields using two separate balanced homodyne detectors, achieving overall collection efficiencies of ηA = 60% and ηB = 77%. This gives a combined measurement efficiency of ηmeas = 58% (see Supplementary Note 5).

By combining slope and dither lock techniques we are able to arbitrarily stabilize the phase of the local oscillators in the range [0, 2π). The photocurrent of each balanced homodyne detector is digitized with a 15 MSa/s analog-to-digital converter. We numerically demodulate the photocurrents at frequency Ω/(2π) = 1.1416 MHz, and low-pass filter the result with bandwidth 200 Hz. (This bandwidth is narrow compared to the mechanical feature of interest, allowing us to apply the infinitely long-filter limit of the model).

Inseparability

We now proceed to characterize the variance of EPR-type variables, as introduced above, to compare with the DGCZ criterion. We choose a common basis θA ≈ θB ≈ Θ ≈ 0 and measure, in sequence, the combinations \(\{{\hat{X}}_{A}^{\Theta },{\hat{X}}_{B}^{\Theta }\}\), \(\{{\hat{Y}}_{A}^{\Theta },{\hat{Y}}_{B}^{\Theta }\}\), and vacuum noise (by blocking the cavity output). Figure 2a, b shows histograms of the measured quadrature data for the X and Y quadratures, along with reference lines for vacuum noise variance in black. Recalling the joint quadrature definitions, we note that the DGCZ criterion involves the diagonal and anti-diagonal variances of the X and Y histograms, respectively. In the figure, we clearly see squeezing in the former, and near-vacuum variance in the latter—thus already indicating violation of the DGCZ criterion. Quantitatively, we find \({\mathcal{I}}=0.83\pm 2 \% ({\rm{stat}}.)\pm 0.3 \% ({\rm{syst}}.)\). The statistical error comes from the number of samples used to estimate the EPR variances and the vacuum noise. The systematic error arises from the estimations of the vacuum noise variance, due to residual classical amplitude noise and mismatch in the photodiode responsivities, at the balanced homodyne detectors (see Supplementary Note 6). We also notice that the variances in the orthogonal directions are at the vacuum level. This does not violate the Heisenberg uncertainty relation, since the pairs of quadratures {\({\hat{X}}_{A},{\hat{X}}_{B}\)} and {\({\hat{Y}}_{A},{\hat{Y}}_{B}\)} commute with each other and are not canonically conjugate observables. We repeat such measurement for different homodyne angles (Θ [−π/2, π/2]) as shown in Fig. 2c. The solid lines are theoretical predictions based on a full model of optomechanical dynamics (taking, in particular, dynamical backaction14 into account, see Supplementary Note 1), using system parameters extracted from fits (see Supplementary Note 5). We find good agreement over all phases, firmly establishing that the two optical modes satisfy the DGCZ inseparability criterion. From Fig. 2c, we also notice that the best two-mode squeezing we achieve is, for the quadrature \({\hat{X}}_{+}\), 1.8 dB below the vacuum noise limit.

Fig. 2: EPR quadrature statistics.
figure 2

a, b 2D histograms of raw X and Y quadrature data, respectively, for Θ ≈ 0. The black dashed circle indicates (2×) the vacuum noise variance, which has a radius of \(1/\sqrt{2}\) (note the axes' scale factor of \(\sqrt{2}\)). The solid black ellipses indicate (2×) the covariance ellipse of the measured data. The black arrows indicate the diagonal/anti-diagonal variances which are relevant for calculating the DGCZ criterion. c Homodyne angle dependence of joint quadrature variances. The purple (orange) dots are the sum (difference) quadrature \({\hat{X}}_{+}\) (\({\hat{Y}}_{-}\)). The average of these yields the DGCZ inseparability criterion (green points). The measurement ensembles contain  ~104 samples, such that the statistical standard error of the variance estimators is  ~2% of the reported values. This contains both the error in the estimation of the EPR variances and the error in the estimations of the vacuum noise variance. The solid line is the theoretical prediction.

Covariance matrix tomography

Having established entanglement, we now quantify it by reconstructing the covariance matrix by Gaussian homodyne tomography. By measuring five different pairs of angles {θAθB} = {0, 0}, {π/2, π/2}, {0, π/2}, {π/2, 0}, {π/4, π/4}, we obtain all necessary intrasystem and intersystem correlations. The reconstructed covariance matrix and theoretical prediction are shown in Fig. 3. From this experimental data, we find a minimum symplectic eigenvalue \(2{\tilde{\nu }}_{-}=0.79\), corresponding to EN = 0.35.

Fig. 3: Covariance matrix of the two optical modes.
figure 3

a Measured (black) and predicted (gray) entries of the covariance matrix. The measurement ensembles contain  ~104 samples, such that the standard error of the variance estimators is  ~2% of the reported values. b Matrix representation of the measured data, to highlight the location of the significant nonzero entries.

Frequency-dependent entanglement

The previous results refer to a sideband quadrature mode at a particular frequency, Ω. We now examine how this entanglement varies as we sweep Ω near the mechanical resonance, Ωm.

(Note that for computational convenience, we do this by calculating noise spectral densities via the Fast Fourier Transform of the raw photocurrents, which corresponds to a mode function \(h(t)=\theta (t)\theta (T-t)/\sqrt{T}\), where θ(t) is the Heaviside function and T ≈ 9 ms is the acquisition time). Figure 4 shows such a frequency-dependent inseparability, as well as its dependence on the homodyne measurement basis, Θ.

Fig. 4: Frequency-dependent entanglement.
figure 4

a Variance of the EPR-type joint quadratures, \({\hat{X}}_{+}^{\Theta }\) (purple) and \({\hat{Y}}_{-}^{\Theta }\) (orange), from sum and difference of the two measured homodyne photocurrents, at angle Θ ≈ 0 and b derived inseparability (green), as a function of center frequency Ω. The solid lines in a and b are fit to a full model (see Supplementary Note 1). The dark gray line represents the minimum symplectic eigenvalue \(2{\tilde{\nu }}_{-}\), obtained by optical homodyne tomography. The modes of the two laser fields are entangled whenever \(2{\tilde{\nu }}_{-} \, < \, 1\). c Theory and d measurements of the inseparability \({\mathcal{I}}(\Theta ,\Omega )\). The green dashed line indicates the measurement shown in b. The horizontal axes are referenced to the bare mechanical frequency Ωm = 2π × 1.139 MHz.

We see that the entanglement criteria can be met for frequencies above and below mechanical resonance, in a manner consistent with the dispersive third term in Eq. (4). The solid lines in Fig. 4a, b are theoretical predictions from the full model, based on a single set of system parameters, obtained from independent measurements or fits to account for drifts (see Supplementary Note 5). Moreover, similar to the measurement in Fig. 3, we can reconstruct the covariance matrix (and corresponding \({\tilde{\nu }}_{-}\)) for each frequency bin, as shown in Fig. 4b. We see that, as expected, \(2{\tilde{\nu }}_{-}\) serves as a lower bound for the inseparability \({\mathcal{I}}\). Since this bipartite Gaussian state is approximately symmetric, from \(2{\tilde{\nu }}_{-}\) we can calculate the entanglement of formation, which is accepted as a proper measure of quantum correlations as a resource15,28,29. Integrating this quantity over a 30 kHz bandwidth yields an entanglement distribution rate of 753 ebits/s.

We emphasize that the optomechanical interaction which generated the entanglement presented above is fundamentally wavelength independent. To illustrate this, we move laser A to ~912 nm, and repeat the measurements of Fig. 4. As shown in Fig. 5, we observe a DGCZ variance below unity and a minimum symplectic eigenvalue \(2{\tilde{\nu }}_{-}=0.92\) for a mode centered at Ω = 2π × 1.142 MHz with bandwidth 915 Hz. The performance is degraded compared to the previous results, due to less efficient light collection at  ~912 nm. Nevertheless, these results establish entanglement of two lasers separated by more than 100 nm in wavelength.

Discussion

In conclusion, we have demonstrated quadrature entanglement of two nondegenerate optical beams via their common radiation–pressure interaction with a mechanical resonator. While applications in optical or microwave quantum communication are conceivable (as realized with entangled light sources based on more traditional optical or microwave parametric oscillators28,30,31,32,33), mechanical platforms offer unique possibilities. In particular, the combination of mechanically mediated microwave15 and optical (this work) entanglement could enable microwave-optical entanglement, based on membrane electro-opto-mechanical systems34,35. This would constitute a much-needed resource for networks of quantum computers based on superconducting qubits. In this context, it is noteworthy that the mechanical interface can in principle also store quantum information. Indeed, for the device presented here, the memory time is ca. 1 ms even for 10 K operation27, easily exceeding storage time in optical fibers.

In our work, entanglement is preserved from the cryogenic mechanical mediator all the way to the laser beams analyzed in room-temperature homodyne detectors. This enhances the prospects of a general class of hybrid quantum systems3,4 based on mechanical interfaces, which could harness entanglement between solid-state (e.g., spin or charge based) quantum systems, typically operating at low temperatures, and itinerant optical fields.

From a more fundamental perspective, it would be interesting to explore concepts at the interface of quantum measurement and entanglement. For instance, the optomechanical interaction in this work can also be interpreted as a strong quantum measurement of the mechanics. This system should be well-suited for studying the usually hard-to-access system-meter entanglement36,37,38.

Fig. 5: Entanglement spanning 120 nm in wavelength.
figure 5

Inseparability \({\mathcal{I}}(\Theta \approx 0,\Omega )\) (green) and minimum of the symplectic eigenvalue \(2{\tilde{\nu }}_{-}\) (gray) from homodyne tomography. The horizontal axis is referenced to the bare mechanical frequency Ωm = 2π × 1.139 MHz.