Stationary optomechanical entanglement between a mechanical oscillator and its measurement apparatus

We provide an argument to infer stationary entanglement between light and a mechanical oscillator based on continuous measurement of light only. We propose an experimentally realisable scheme involving an optomechanical cavity driven by a resonant, continuous-wave field operating in the non-sideband-resolved regime. This corresponds to the conventional configuration of an optomechanical position or force sensor. We show analytically that entanglement between the mechanical oscillator and the output field of the optomechanical cavity can be inferred from the measurement of squeezing in (generalized) Einstein--Podolski--Rosen quadratures of suitable temporal modes of the stationary light field. Squeezing can reach levels of up to 50% of noise reduction below shot noise in the limit of large quantum cooperativity, corresponding to a back-action limited position sensor. Remarkably, entanglement persists even in the opposite limit of small cooperativity. Viewing the optomechanical device as a position sensor, entanglement between mechanics and light is an instance of object--apparatus entanglement predicted by quantum measurement theory.

It has been a long-standing prediction by Genes et al. [16] that optomechanical entanglement persists under certain conditions in steady state under continuous-wave drive, see Ref. [17] for a comprehensive review.As explained in Ref. [18], the experimental verification of this effect is very challenging, because the state of the mechanical oscillator is not accessible directly.Measuring limit, where the linewidth of the optomechanical cavity is large compared to the mechanical frequency.We consider this scenario because it allows for a comprehensive and insightful analytical characterization of the problem.We will show that entanglement can be detected in the form of squeezing of Einstein-Podolski-Rosen modes of light, or suitable generalizations thereof.We predict squeezing to scale inversely proportional to the quantum cooperativity and asymptotically reach 50% of noise reduction below the shot noise level.This leaves a safe margin for experimental imperfections, such as finite detection efficiency.Remarkably, we find that optomechanical entanglement persists even for a quantum cooperativity below unity, albeit only to a small extent.
The operating regime of the scheme we propose corresponds to the generic configuration of an optomechanical position or force sensor.In a continuous position sensor the measured object (the mechanical oscillator) and the measurement apparatus (the light field) share entanglement.This is the central idea of the quantum mechanical measurement theory [24][25][26]: the process underlying a physical measurement is ultimately an entangling interaction between object and the measuring apparatus (observer).This applies regardless of whether the measurement is performed at high quantum cooperativity, for which it is limited by measurement back-action, or at low cooperativity, where it is limited by shot noise or thermal noise.Under normal conditions, it is virtually impossible to detect the entanglement between object and apparatus as it usually involves an uncontrollable variety of environmental degrees of freedom due to the amplification associated with the measurement.It is intriguing to see that with an optical-mechanical sensor it is possible to capture this entanglement in a feasible measurement and thus shift the Heisenberg-von Neumann cut between object and observer.
This article is organized as follows: in Sec.II we introduce necessary elements of optomechanics.Our argument and scheme to detect optomechanical entanglement is developed in Sec.III.We present our results and predictions regarding detectable signatures of optomechanical entanglement in Sec.IV.

II. ELEMENTS OF OPTOMECHANICS
We consider a standard optomechanical system [27] comprising a single mechanical (oscillatory) mode interacting with a single light mode of a cavity.Mechanical and light fields are bosonic, described by two pairs of dimensionless Hermitian operators, x m /p m and x c /p c , referred to as position/momentum and amplitude/phase quadratures of the mechanical and cavity mode, respec-tively.These operators obey equal-time canonical commutation relations, where α, β = m, c and δ m,c is the Kronecker delta symbol.Throughout this work = 1.We will also use ladder operators c α and c † α defined by obeying [c α (t), c † β (t)] = δ α,β .In a frame rotating at the frequency ω d of the field driving the optomechanical cavity, the Hamiltonian of the linearized dynamics is given by [27] We denote the oscillation frequency of the mechanics by ω m , ∆ ≈ ω d − ω c is the detuning of the drive field from cavity resonance, and g is the coupling strength (depending on the drive power).The interaction contains two processes that create sidebands in the cavity field: the down-conversion of a phonon and a photon at a lower energy than the drive (that is, Stokes scattering to sideband frequency −ω m ), and the state-swap of a phonon onto a photon at higher energy than the drive (that is, anti-Stokes scattering to sideband frequency +ω m ).The cavity mode is coupled to the free electromagnetic field outside the cavity.The ladder operators describing this one-dimensional field, c in (t) and c † in (t), have the following Markovian correlation functions We neglect here thermal occupation numbers for optical frequencies.
The thermal bath for the mechanical oscillator is modeled by quantum Brownian motion damping with the Hermitian noise operator ξ [28,29].At high temperatures of the mechanical bath n th ≈ k B T bath / ω m 1 and/or large mechanical quality factor Q = ω m /γ m 1 its correlators are approximately Markovian, reflected by Including cavity (power) decay at rate κ and viscous damping of the oscillator at rate γ m , the open-system dynamics is described by the quantum Langevin equations (QLE) [27] [30] ẋm = ω m p m , (6a) In the following we consider the special case of a resonantly driven cavity, ∆ = 0, which corresponds to the standard configuration of an optomechanical position or force sensor.The generalisation to nonzero detuning is straightforward, but results in significantly more involved analytical expressions that we will not reproduce here.Note that x in and p in in Eqs.(6) correspond to shot noise only, because we work in a suitably displaced frame where mean amplitudes were shifted to zero, see Ref. [27] for details.
Our study always assumes stable dynamics such that the steady-state is reached eventually [31].Equations ( 6) can be solved easily in Fourier space [16]; see App.B 1 for the Fourier transform conventions we have adopted.The quadratures of the light emitted by the cavity are given by input-output relations [32] and take the form We have introduced here the mechanical and the optical susceptibilities and the reflection phase We assume Q = ω m /γ m 1 as a typical feature in micro-optomechanical setups.This allows to approximate the poles of the mechanical susceptibility, Eq. ( 8), by ω ± ≈ ±ω m − iγ m /2, such that Both are strongly peaked at ±ω m .Close to these frequencies and in the bad cavity regime (κ ω m ), which will become relevant later, we approximate Eqs. ( 9) and (10) by, respectively, In these approximations, the characteristic response time 1/κ of the intra-cavity field is the shortest time scale of the system, such that the intra-cavity field is adiabatically eliminated from the dynamics.In other words, the mechanical oscillator is effectively directly coupled to the output field, without any spectral filtering due to the cavity, cf.Eq. (11c).In this limit, one can rewrite Eq. (7b) as where we have defined the readout rate The first term of p out in Eq. ( 12) is the contribution of shot noise reflected by the cavity.The second term is the contribution from shot noise driving the mechanical motion and being transferred back to the light via the interaction -this effect is called back-action.The last term is the contribution of the thermal fluctuations from the mechanical bath, mapped onto the light via the optomechanical interaction.In the context of position or force sensing, this term constitutes the signal to be determined from a measurement of the phase quadrature p out .The relative size of the readout term to the thermal noise in Eq. ( 12) is the so-called quantum cooperativity The last approximation holds in the high temperature limit, and we defined the thermal decoherence rate of the mechanical oscillator which sets the scale of the thermal force in Eq. (6b).At a fundamental level, the transduction of information on position or force is associated with the generation of correlations, or even quantum entanglement, between the observed object (here, the mechanical oscillator) and the measurement apparatus (the light field).In the next section we will explain how stationary optomechanical entanglement can be proven unambiguously based solely on measurements of the light field.Consider two consecutive pulses of light, an early pulse (E) and a late pulse (L), whose temporal envelopes have support on non-overlapping time intervals [t E − τ, t E ] and [t L , t L + τ ], respectively, with t E < t L and pulse duration τ > 0. Assume that the two pulses and the mechanical oscillator (M) are initially in a product state ρ init EML .Furthermore, let the dynamics be such that the mechanical oscillator interacts first with the early pulse according to a quantum channel (i.e., a completely positive map) E EM , followed by an interaction of the oscillator with the late pulse described by E ML , see Fig. 1.The final state of the three systems is ρ fin EML = E ML E EM (ρ init EML ) .One sees that if the reduced state for the mechanical oscillator and the early pulse ρ EM = tr L (ρ fin EML ) is separable, then the state of early and late pulse ρ EL = tr M (ρ fin EML ) is separable too.Conversely, if ρ EL is entangled, ρ ME must have been entangled.Thus, the measurement of entanglement between early and late pulses implies optomechanical entanglement.A detailed presentation of this argument is given in App. A.
Based on this reasoning, entanglement between a mechanical resonator and a travelling wave pulse of a microwave field had been demonstrated in Ref. [14], along the lines outlined in Ref. [18].Optomechanical entan-glement had been deduced from the measurement of entanglement between two consecutive microwave pulses that originated from sequential interactions with the mechanical oscillator.In the scheme of Refs.[14,18], the two processes E EM and E ML are different in nature, and constructed to optimize entanglement generation in E EM and state readout during E ML .This requires, in particular, to work in the resolved-sideband limit (κ/ω m < 1), and to set the driving field blue-detuned (∆ = ω m ) in the first (E-M) interaction and red-detuned (∆ = −ω m ) during the second (M-L) interaction.In both processes the applied fields are pulsed and steer the optomechanical system through non-stationary, non-equilibrium states.
In the following we will show that the same logic for verifying optomechanical entanglement can be applied in a regime in which the optomechanical system is driven by continuous wave light, the detuning is fixed, and the mechanical oscillator and the output field of the cavity are in a stationary state.The two successive temporal pulses, which are crucial for our argument, are only extracted in post-processing from the continuous measurement on the output field.Because the overall state of the mechanical oscillator and the output field is stationary, it does not matter which intervals [t E − τ, t E ] and [t L , t L + τ ] are extracted from the stationary homodyne current: any of its properties, such as entanglement of early and late pulses, will depend only on the pulse length τ and separation T sep = t L − t E between the pulses, but not on the particular instants of time t E , t L .Formally we define annihilation operators corresponding to these temporally-ordered modes by where f E (t) and f L (t) are the temporal mode functions of the early and late modes, respectively.The properties of c out are determined by Eqs.(7), which describe the stationary state of light in conjunction with the properties of the ongoing noise processes Eqs. ( 4) and (5).
Because the two modes E and L are ordered in time and interact sequentially with the mechanical oscillator, the scenario of Fig. 1 applies.Thus, optomechanical entanglement is demonstrated trough detecting the two modes and showing that they share entanglement.
It is instructive to compare this once again to the pulsed scheme of Refs.[14,18].There the properties of early and late pulses had to be inferred from an integration in time of the respective equations of motion, Eqs.(6), with different detuning ∆ = ±ω m for E and L.Here we can infer all properties of E and L from the stationary solutions, Eqs.(7), for a fixed detuning.In the pulsed scheme the description of the three modes E, L and M was essentially complete in the sense that the dynamics were designed so that no correlations to any other mode of the light field were established.Here, the restriction to the two pulses E and L is a massive simplification: in reality, the mechanical oscillator will exhibit correlations to more modes than E and L, and the stationary light field will exhibit a large variety of internal correlations, e.g., in the form of ponderomotive squeezing.Still, for demonstrating stationary entanglement of the mechanical oscillator and light, it will be sufficient to consider the three modes E, M and L.
For the following discussion it will be convenient to formally allow for infinitely extended pulses, corresponding to τ → ∞, whose pulse envelopes f E/L (t) tend to zero for t → ∓∞, respectively.The characteristic scale at which the envelopes tend to zero defines an effective pulse length.Furthermore, because of stationarity, we can set arbitrarily t E to −T sep /2 and t L to T sep /2 without loss of generality.We write the mode operators r i , for i = E, L, as where f E (t) has to be an anti-causal function with respect to −T sep /2 (i.e., f E (t) = 0 for t > −T sep /2) and f L (t) has to be causal with respect to T sep /2 (i.e., f L (t) = 0 for t < T sep /2).The r i must fulfil bosonic commutation relations, [r i (t), r † j (t)] = δ i,j , which correspond to an orthonormalisation constraint for the mode functions With Eq. ( 2) we obtain the quadrature operators x i , p i of the early and late modes of the output field.
In order to prove stationary optomechanical entanglement, we need to identify suitable mode functions and system parameters that result in an entangled state ρ EL of early and late pulses, which can be verified from a measurement of the output light field of the optomechanical cavity.
In the following Sec.III B we motivate a choice of entanglement criterion.This specifies what information the measurement procedure must retrieve from the output field.In Sec.III C we motivate a particular class of temporal mode functions to process measurement data.This choice is based on our understanding of the stationary optomechanical dynamics and leads to choosing the particular parameter regime of unresolved sidebands and undetuned drive.These choices determine an implementable experimental scheme that we can study analytically; the results of the study are presented in Sec.IV.

B. Entanglement criteria
We note first that the optomechanical dynamics is linear and the driving fields correspond to Gaussian white noise.This implies that the stationary state of the optomechanical system and the output field is Gaussian, and that the state of the two temporal modes in the output field, introduced in Eq. ( 17), is Gaussian too.First moments and covariances of canonical operators x i , p i (i = E, L) determine Gaussian states completely.As a consequence all entanglement properties are fully determined by the latter.Entanglement of the two-mode Gaussian state of E and L can be quantified by means of suitable entanglement measures [33][34][35] (such as the Gaussian entanglement of formation [36] or the logarithmic negativity [37][38][39]), which can be calculated from the experimentally accessible covariances of the state ρ EL .What is more, such entanglement measures often provide bounds to non-Gaussian states for given covariances.In order to proceed in a fully systematic way, one could try to optimize temporal mode functions in Eq. ( 17) with respect to an entanglement measure and other system parameters (such as optomechanical coupling g, etc.).However, due to the highly nonlinear dependence of entanglement measures on the state ρ EL such an optimization seems unfeasible.
Instead, we make use of entanglement criteria that are linear in covariances, so second moments of quadratures [40][41][42].The geometry underlying these criteria is that they specify hyperplanes in the set of covariances: they are separating hyperplanes from the convex set of covariances that are compatible with separable Gaussian states [42].That is to say, on the level of covariances, the geometry is basically that of an entanglement witness for quantum states.Since each test is linear in the covariances, each test defines an implementable and feasible measurement procedure.What is more, one can argue that a strong violation of such criteria is accompanied by a quantitative statement, in that several entanglement measures are lower bounded by quantitative violations of such entanglement tests [42][43][44][45].Moreover, for an anticipated covariance, one can efficiently find the optimal test that best certifies the entanglement present in a Gaussian state [42].
For our analytical consideration we will resort to a special case of such a criterion, referred to as Duan's criterion [40]: the so-called EPR-variance (∆ EPR ) where Our notation for the variances is where the expectation values are taken with respect to the initial density operator It quantifies the simultaneous correlations between pairs of quadratures of different modes.For two-mode Gaussian states, implies (only sufficient) that the state is entangled.In terms of ladder operators, the EPR-variance can be written as This relation can be readily evaluated from a record of homodyne measurements, along with appropriate postprocessing in order to extract appropriate (anti-)causal temporal modes.Again, the Duan criterion is in general not an optimal test, but can be optimized for a given Gaussian state.

C. Temporal modes
In this section, we address the question of how to choose the mode-functions f E/L (ω) in Eq. ( 17) in order to achieve largest violation of the separability bound Eq. ( 21).The EPR-variance can be expressed as a quadratic form of the mode-functions, which we show explicitly in App.B, cf.Eq. (B8).Unfortunately the minimization does not map to a simple eigenvalue problem due to the required (anti-)causality and normalization of the mode functions.
We therefore proceed by choosing a suitable variational class of mode functions that we motivate as follows: in section II, we explained that the optomechanical unitary dynamics produce sidebands around the driving frequency at ±ω m (in the frame rotating at the drive's frequency).The lower (red) sideband is produced by the down-conversion part (also termed two-modesqueezing) of the interaction and it is known to produce EPR-like entanglement [16,18]: this sideband is entangled with the mechanics.The upper (blue) sideband is produced by the beam-splitter contribution of the interaction and corresponds to a coherent state-swap of the mechanical and cavity fields (preserving quantum correlations with the red sideband): this sideband encodes the state of the mechanics, which was entangled with the red sideband some moments before.The blue sideband is thus entangled with the light in the red sideband at previous times.If ∆ = 0 and if the cavity linewidth κ is not much narrower than ω m , then both sidebands will have significant spectral components in the measured outputthis is the reason for choosing the regime of unresolved sidebands and resonant drive.We choose early and late mode functions that capture the temporal order of these processes: the early mode extracts information from the red sideband and the late mode from the blue sideband.This is done with a demodulation of the detected signal by e ±iωmt , which makes the corresponding spectral component resonant.
We expect the quantum correlations at different points in time to decrease with the duration separating them.This is because the incoherent coupling to the mechanical bath eventually destroys coherent correlations between mechanics and light.We therefore opt for pulses with an effective duration that is smaller than the thermal decoherence time Γ −1 th , cf.Eq. ( 15).As a variational class of mode functions we choose exponentially decaying envelopes for the demodulating phases, see Fig. 2, where θ is the Heaviside step function and N = (2Γe ΓTsep ) 1/2 is the normalization constant.These mode functions are anti-causal and causal, respectively, and they are orthonormal by construction.The pulse decay rate Γ is a variational parameter with respect to which we will minimize ∆ EPR .

D. Experimental task
Sec. III defines a readily implementable scheme to detect optomechanical entanglement between temporally- ordered modes of the output field.The optomechanical cavity is driven on resonance in the non-sidebandresolved regime by a constant input laser while the cavity output field is continuously monitored.The choice of the exact measurement scheme depends on the entanglement criterion of interest.Here Eq. (19) [or equivalently Eq. (22)] prescribes what information the detection scheme in an experiment must provide.For exam-ple, the equal-time correlators in Eq. ( 22) can be evaluated from records of two homodyne measurements of the phase and the amplitude of the output field.From such a record, temporally-ordered modes can be extracted according to Eq. ( 17).The detection of entanglement between the temporally-ordered modes implies that the mechanical oscillator and the output light share stationary entanglement.Figure 3 summarizes this procedure.

IV. RESULTS
This section has two parts: first we provide an approximate formula for the EPR-variance and comment on what it teaches us about entanglement in our system.In the second part we benchmark the approximate formula with its exact counterpart.The latter is very cumbersome in its symbolic form and on its own provides only numerical results, whereas our approximate formula gives generic behaviour with respect to the parameters.

A. Approximate formula for the EPR-variance
Using the mode functions of Eqs.(23) in the definition of the mode operators in Eq. ( 17), the EPR-variance in Eq. ( 22) can be evaluated by means of the input-output relations Eqs.(7) and the noise correlators Eqs. ( 5) and (4).A detailed derivation is given in App.B.
Here we give the final result, which is compared to (and backed up by) exact numerical calculations in Sec.IV B, The readout and thermal decoherence rate, Γ ro and Γ th , have been defined in Eqs. ( 13) and (15).We included a finite detection efficiency η ≤ 1. Equation ( 24) is derived for a resonant drive ∆ = 0 and the regime κ ω m γ m such that approximations Eqs. ( 11) apply.Furthermore, in order to arrive at an analytically tractable formula, we made the technical assumption that the pulse envelopes decay at a rate Γ fulfilling Here C q is the optomechanical quantum cooperativity, cf.Eq. ( 14), and C cl := 4g 2 /κγ m is the classical cooperativity.The approximation assumes the large temper-ature limit.Equation (25) has to be read as a limitation on the optomechanical coupling g.As we will see, for couplings that violate Eq. ( 25) the EPR-variance is not a suitable entanglement criterion anymore.Thus, in the regime in which the EPR-variance is of interest, Eq. ( 25) is a self-consistent restriction.
It is useful to characterize the circumstances for which ∆ EPR is minimal, so that the violation of the separability bound, Eq. ( 21), is least susceptible to inevitable experimental imperfections.To this end, we set φ = 0 in Eq. (24).The first term in the square brackets is always positive and monotonically decreasing with the decay rate of the temporal modes Γ.The second term is negative and monotonically increasing with Γ.Only if the second term can compensate the first one is it possible to detect entanglement.In view of this and the overall scaling in front of the square brackets with respect to Γ, it is clear that there is an optimal pulse bandwidth yielding an optimal EPR-variance.
We discuss first the result of this minimization for ideal detection, η = 1, and vanishing pulse separation, that is, in zeroth order of γ m T sep .A straightforward calculation gives for which the approximation holds when n th 1.We expected that the pulses have to decay faster than the thermal decoherence rate in order to preserve the coherence between the pulses with large probability.This means that the pulses' bandwidth can not be too narrow.The factor of 2 accounts for the fact that the coherence must be kept between both pulses, i.e., there may be no significant incoherent leakage over times ∼ 2/Γ opt .Increasing the readout rate (or the cooperativity) leads to shorter optimal pulses.We interpret this tendency as a tentative to extract information from the output light as fast as the dynamics allows, in order to minimize the chances of decoherence.This seems a natural strategy because, in our model, thermal decoherence is the only channel through which entanglement can be lost.However, this might not be the best strategy in a real experiment: as pulses become shorter, their bandwidth increases and spectral incoherent features of technical noise and/or additional mechanical modes will be resolved by the pulses.We expect that this will prevent observation of entanglement if not accounted and controlled carefully.Moreover, the form of the optimal bandwidth Γ opt is only valid when Eq. ( 25) holds, which limits the maximal bandwidths we can predict this way.The optimal pulse bandwidth in Eq. ( 26) is compatible with the assumption in Eq. ( 25), if 2(C cl +n th ) 3/2 Q.For large optomechanical coupling g, which violate this condition, entanglement will not be detectable in the form of EPR-squeezing and one should instead consider more general entanglement witnesses, as we will see.
Formula (24), for the optimal choice of bandwidth Γ opt , predicts Surprisingly, entanglement persists for arbitrary small but nonzero cooperativity.We remind the reader that the EPR-variance refers to directly measurable entanglement between temporally-ordered modes of light, which witnesses stationary entanglement between the mechanical oscillator and the output field of the optomechanical cavity.Equation ( 27) constitutes the main result of our work.For large quantum cooperativity ∆ EPR approaches 50% of squeezing below shot noise.This limit hints at the fact that we consider only two (out of infinitely many) light modes that encode correlations with the mechanics.When finite detection efficiency and pulse separation are taken into account, the minimal EPR-variance, Eq. ( 27), is modified as follows (to first order in γ m T sep ) In the relevant regime of large cooperativity, C q > 1, we have Γ ro > Γ th such that the pulse separation T sep needs to stay well below the thermal decoherence time Γ −1 th , as one would expect.
From the derivation detailed in App.B, it is possible to give a detailed account on the physical origin of each term in Eq. ( 24).The factor 2 at the front is the shot noise contribution of the light field that we would observe in the absence of optomechanical coupling (i.e., if Γ ro = 0).The first term in the square brackets is due to auto-correlations in thermal noise and in back-action noise, contributing to the two pulses.This term has two parts: a non-negative part coming from intra-pulse correlation in each of the two temporal modes (early or late with themselves), and a negative part coming from the inter-pulse correlation between early and late mode.The net effect of both parts is always positive, therefore auto-correlations in thermal noise and back-action noise do not contribute to entanglement.The last term in the square brackets corresponds to cross-correlations of shot noise and back-action noise.These correlations are also at the basis of ponderomotive squeezing typically observed in the frequency domain [4,5].Here, the correlations refer purely to inter-pulse correlation of the two temporal modes; the intra-pulse correlation gives an exact zero contribution, as detailed in App.B. In this sense, an EPR-variance below 2 is due to ponderomotive squeezing between the early and late modes.

B. Comparison with numerical results
Ref. [16] provides a procedure to express analytically (in integral form) the covariance matrix of arbitrary temporal modes of the continuous steady-state output field of an optomechanical cavity.Arbitrary parameter regimes (e.g., nonzero detuning, good cavity limit, strong coupling, etc.) can be studied numerically this way, as long as the system is stable and reaches a steady-state.The form of the mode functions defined in Eqs.(23), for van-  27) (red solid line), exact EPR-variance (green circles) and optimal entanglement witness (blue dots) versus optomechanical coupling g (lower x-axis) and quantum cooperativity Cq (upper x-axis).Threshold for separable states is at 2 (solid black line), the yellow shaded area represents the region of entanglement for the EPR-variance and the optimal witness.This applies to all other figures.
ishing T sep , allows one to compute symbolically the exact expression of the covariance matrix, and the EPRvariance is readily obtained from the entries of the latter.The exact formula for the EPR-variance is cumbersome and not very informative in its symbolic form, but it provides numerical results upon evaluation with fixed parameters.We compare below the approximate formula, Eq. ( 24), to this approximation-free method.
In the following, unless otherwise stated, we assume in all figures ω m /2π = 1 MHz, κ = 10 ω m , Q = ω m /γ m = 10 8 , and n th = 10 4 .This parameter set is partly inspired by a recent experiment [46].All plots refer to the case η = 1 and T sep = 0 (which corresponds to the limit where T sep does not significantly alter ∆ EPR ).
Figure 4 shows the minimal EPR-variance versus optomechanical coupling g (and quantum cooperativity C q ), comparing the approximate formula, Eq. ( 27) (red line), to the exact results (green circles) of the approximation-free symbolic expression.In the latter, for each g, we sweep over Γ to determine the minimal value of ∆ EPR .For cooperativities C q 3 the plots overlap well.In particular, the exact result confirms the surprising presence of entanglement at small cooperativities.A discrepancy between the approximate formula and the exact result emerges because the restriction on the pulse bandwidths, Eq. ( 25) for Γ = Γ opt , does not hold and the formula becomes inaccurate.The exact ∆ EPR (green circles) rises above the separability threshold 2 when the cooperativity increases: this means that,   27) (solid lines), exact EPR variance (circles), optimal entanglement witness (dots) versus pulse bandwidth Γ for cooperativities of Cq = 0.1, 1, 10 in red, green, blue, respectively.
in the large cooperativity regime, the entanglement is not in the EPR-mode anymore.For all two-modes Gaussian states it is possible to find the optimal entanglement witness, based on linear combinations of second moments, which decides whether the state is entangled [42][47].
Because the EPR-variance is a (in general sub-optimal) witness based on second moments, we can match the separability threshold for any witness to the one of the EPR-variance and plot the optimal witness on the same y-axis.The blue dots in Fig. 4 are the optimal-witness values and they monotonically decrease as C q increases, hence, in principle, increasing C q does not make the detection of entanglement more difficult.This is an example of ∆ EPR being sub-optimal, to the point that it cannot detect entanglement.Interestingly, the approximate formula for the EPR-variance and the optimal witness overlap well.We interpret this to be a result of the general scaling of squeezing as 1/C q in the limit of large cooperativity.
Figure 5 shows how the EPR-variance changes with the pulse bandwidths Γ according to the approximate formula, Eq. ( 24) (solid lines), and compared to the exact results (circles).All curves display a single global minimum that is not a sharp feature.Thus the optimisation over the pulse bandwidths is not too difficult in principle; especially in an experimental scenario where all the parameters are not known exactly.According to formula (24), larger pulse bandwidths will always yield some entanglement though one has to keep in mind the restriction of Eq. ( 25) that limits the validity of the formula with respect to large Γ.As mentioned in the previous section, thermal noise on the mechanics is the sole decoherence channel of our model, therefore pulses decaying faster  27) (red solid line), exact EPR variance (green circles), optimal entanglement witness (blue dots) versus mechanical bath temperature, parametrized in terms of the mean thermal occupation number, for fixed optomechanical coupling g = 2π × 15.8KHz.
than the thermal decoherence time 1/Γ th will display some coherence.In practice, broadband pulses (short in time) will resolve other incoherent spectral components (electronic filters, additional mechanical modes, etc.) not accounted for in the present model.Close to the minima, we see good agreement for C q 1, consistently with Fig. 4. The exact EPR-variance attains larger values at larger cooperativity because it is a sub-optimal witness, as in Fig. 4. On the curve C q = 1, the approximate formula and the exact results agree across the minimum and no more at larger bandwidths where the restriction of Eq. ( 25) does not hold.
In Fig. 6 we illustrate the dependence of the minimal EPR-variance and the optimal witness on the bath temperature of the mechanical oscillator, parametrized by its mean thermal occupation n th .The cooperativity increases like the inverse of n th .The approximate formula, Eq. ( 27) (red line), and the exact numerical result (green circles) agree when C q > 1.The exact ∆ EPR and the optimized witness appears to saturate at a minimal value of 50% of squeezing below shot noise, just like the approximate formula predicts.Formula (24) becomes inaccurate as n th grows large where, again, the restriction of Eq. ( 25) does not hold any more.At large temperature of the mechanical bath, the optimal witness (blue dots) predicts retrievable entanglement, which is similar to the behavior at small coupling on Fig. 4.
In realistic experiments, zero detuning is a challenging regime of operation because it is at the border to the unstable regime of the optomechanical dynamics at blue detunings.Therefore, in practice, the laser drive is usu- .Exact EPR variance (circles) and optimal entanglement witness (dots) versus red detuning of the driving field from cavity resonance.The optomechanical coupling is kept constant corresponding to a quantum cooperativity of Cq = 0.1, 1, 10 for red, green, and blue symbols, respectively.The dashed lines are the approximate minimal EPR-variance formula, Eq. ( 27), valid only at ∆ = 0 (displayed for comparision).
ally slightly red detuned from the cavity resonance.In Fig. 7 we study the effects of detuning on the minimal EPR-variance with the exact expressions of the covariance matrix: the circles correspond to the exact ∆ EPR and the dots to the optimal-witness values.For comparison, the values of formula Eq. ( 27) (only valid for zero detuning) are displayed as dashed lines.As the detuning approaches the cavity linewidth κ the separability bound violation diminishes because the strength of the sidebands' spectral components in the cavity output field is reduced by the cavity profile.This is analogous to decreasing the detection efficiency (passive losses), in the sense that the spectral weight of the signal extracted by the mode functions decreases relative to the shot noise component.On the other hand, detuning does not introduce any additional noise, therefore the EPRvariance does not increase above 2.In Figs. 4 and 5, where the detuning was zero, the exact EPR-variance increases with cooperativity and Fig. 7 shows this too for large cooperativities and small detunings.Remarkably, at large cooperativity, some detuning reduces the minimal EPR-variance down to an asymptotic value close to the prediction of the approximate formula Eq. ( 27) for zero detuning.

C. Remarks on shot-noise levels and multi-mode entanglement
In this section, we comment on the -somewhat subtle -role of shot noise levels in experimental entanglement detection in optomechanical systems.Given the canonical coordinates O = (x L , p L , x E , p E ), one can capture the second moments in a 4 × 4 covariance matrix Ξ with entries for vanishing means or first moments O j for j = 1, . . ., 4. Any covariance matrix of a quantum state satisfies the Heisenberg uncertainty principle [33][34][35] Ξ + iσ ≥ 0, where σ is the symplectic matrix incorporating the canonical commutation relations.
We will now turn to entanglement criteria.Any entanglement test linear in Ξ can be written in the form where X is the 4 × 4 matrix that captures the measurement settings.It takes a moment of thought -and is explained in Ref. [42] -that the Duan criterion of Ref. [40] corresponds to a specific such choice for X.Such tests are not only convenient for their simplicity in experimental implementations: their linear nature in Ξ renders the assessment of confidence regions of tests simple and feasible.
There is a subtle aspect that one needs to keep in mind when experimentally verifying entanglement: for a bipartite Gaussian state, a state being entangled and having a positive partial transpose is one-to-one.That is to say, exactly the covariance matrices of separable Gaussian quantum states will correspond to covariance matrices that satisfy where Ξ Γ is the covariance matrix of the partial transpose.A slightly inaccurate experimental assessment of the shot noise level will correspond to a matrix xσ where the real x is slightly different from unity.That means that when in such a picture the matrix Ξ is slightly unphysical, it seemingly comes along with the Gaussian state being certified as being entangled even if it is not.Therefore, to certify entanglement with confidence it is crucial to obtain a precise shot noise characterisation, taking into account statistical fluctuations and any temporal drift.Note that these effects can be highly relevant even if the states involved are close to pure.In yet other words, if covariance matrices recovered are close to being pure [and hence close to the boundary of the convex set of covariance matrices defined by Eq. ( 33)], enormous care is necessary when drawing the conclusion that the state is entangled: this may be an artefact of a not sufficiently determined shot noise level.This issue is aggravated in the multi-mode case, which is specifically interesting when assessing entanglement between several optical and mechanical modes.Entanglement tests still take the form of Eq. ( 31), now Ξ being a symmetric 2n × 2n matrix for n modes in total [42].Such tests are highly convenient for detecting multi-mode entanglement in optomechanical systems, in that the optimal test for an anticipated quantum state can be efficiently found by means of methods of semidefinite programming [42], so convex optimization techniques [48].
The above mentioned situation is now expected to be generic: in common experiments, some modes of local subsystems will be close to being pure.A more sophisticated statement of this kind is that in similar ways as natural quantum states are close to being low-rank states [49], commonly encountered covariance matrices Ξ will have a low quantity r that can be seen as a "symplectic rank": if m is the number of unit eigenvalues of −(Ξσ) 2 , then A pure Gaussian state will have r = 0, a state each mode of which is mixed r = n.In this language, commonly encountered Gaussian states will feature covariance matrices that are well approximated by covariance matrices that have a value of r different from n: so some modes will be very susceptible to inaccurately estimated shot noise levels.If they are slightly underestimated, then most naturally encountered states will seem entangled, even if there is no entanglement in the system.This is an important observation little appreciated so far in the literature.This observation also comes hand in hand with the fact that in practical recoveries, a slight unphysicality of the covariance matrices goes hand in hand with entanglement being detected.Having said that, this applies only to the situation of inaccurately assessed shot noise levels: in some setting, even a weak signal in optomechanical systems subjected to very high noise levels can lead to strongly significant predictions [50].Again, this is an aspect one has to be very careful about when assessing optomechanical entanglement.

V. CONCLUSION
In this work, we have presented an argument to demonstrate the presence of stationary optomechanical entanglement in a continuously driven optomechanical cavity.The argument relies on the extraction of temporally-ordered modes from the continuous photocurrent of homodyne detection in post-processing -no modulation of the drive is needed.The analytic study of a specific scheme shows that EPR-squeezing can reach up to 50% below the threshold of separable states for large quantum cooperativity C q .Our approximate formula for the EPR-variance predicts that this limit is approached as C −1 q .Remarkably, this limit seems to hold also for the exact EPR-variance as well as for optimal entanglement witnesses which are linear in the second moments.It would be interesting to establish this as a strict limit.
We have studied a specific class of variational mode profiles, which are physically well motivated, and are theoretically capable to reveal entanglement.Nevertheless, a systematic optimization of the mode functions would be desirable.This could be based, e.g., on the approach indicated in the end of App.B 2.
The single-mode model of the mechanical oscillator has been sufficient to fully develop our entanglement verification argument.Now, single-mode resonators are rare and it turns out that our scheme can be generalized to multiple mechanical modes: for each modes an early and late pair of mode operators, Eq. ( 17), are defined such that the covariance matrix, of dimension twice the number of modes included, can be computed.Entanglement in this state can be assessed with respect to the bipartition early-late.The inclusion of more modes improves the purity of the reconstructed state, in the sense that spectral contributions from the included modes are treated as signal rather than strong noise.This procedure might significantly improve experimental attempts.The original and detailed study has been presented in Ref. [51].
We hope that our scheme provides a viable pathway to demonstrate stationary optomechanical entanglement, as predicted over a decade ago [16].We think that the demonstration of this effect would constitute a textbook experiment exploring the Heisenberg-von Neumann cut between measured system and measurement apparatus in the example of an optomechanical position meter.

A. APPENDIX: LIGHT-LIGHT ENTANGLEMENT IMPLIES LIGHT-MECHANICS ENTANGLEMENT
In this Appendix we provide the details of the indirect proof of entanglement in the scenario shown in Fig. 1.The initial state is a product state of E, M, and L, The first quantum channel E EM acts non-trivially on E and M only.If this quantum channel does not generate entanglement between E and M, then In this case, any subsequent channel acting on M and L cannot generate entanglement between E and L: the final state is and the reduced state of E and L is where ρi L = tr M E ML (ρ i M ⊗ ρ L ) .The state ρ fin EL is separable.In essence, this statement is equivalent to the obvious fact that entanglement cannot be generated by local operations.In section III A we use the contraposed statement: if ρ fin EL is entangled, then E EM must have created E-M entanglement.
In the literature we find that it is possible to distribute entanglement between parties E and L with a mediator M, in such a way that M is never entangled with E and L [19][20][21][22][23].This is not a contradiction with our statement because we assume that the initial state is uncorrelated.shorthand notations and symmetries of Eqs.(B4).The thermal noise integrals need to be combined appropriately to take the from of Eqs.(5) (in frequency space).Each correlator at this point is a single frequency integral and the EPR-variance, Eq. ( 22), can be written as the following explicitly real expression This relation is exact if the drive is not detuned with respect to the cavity.It also reveals how the driving noises couple to each other and to the mode functions: S(ω)S(−ω), B(ω)B(−ω) and T (ω)T (−ω) come from the auto-correlations of the noises, while S(−ω)B(ω) comes from the correlation between shot noise and backaction noise.The correlators are overlapped either with f j f * j -we call them intra-mode integrals -or with f j f k -which we call inter-mode integrals.
So far, we did not use that f E and f L are time-ordered and non-overlapping, therefore Eq. (B8) is valid for arbitrary mode functions satisfying the orthonormality constraint, Eq. (18).Moreover, Eq. (B8) takes a compact matrix form.Indeed, notice that S(ω)S(−ω) = 1, and call It is interesting to look at the problem of minimising Eq. (B10) under the constraints that the mode functions are time-ordered -necessary to infer optomechanical entanglement based on light modes entanglement -and orthonormal, Eq. ( 18).Titchmarsh's theorem [52,53] states that the causality constraint is equivalent to mode functions determined by their real parts only.Therefore, the constraints can be formulated in a quadratic form such that the quadratic minimisation problem can be solved in polynomial times with linear programming methods, according to [54].We were not successful with a simple and naive approach to perform the minimisation this way.We attempted to solve it as an eigenvalue problem but were not successful either: the main issues were the wide parameter spread (ratio between largest and smallest scale is 10 8 ) and the form of the normalisation constraint.Moreover, optimising the mode functions to minimize ∆ EPR is sub-optimal but we expect that repeating our analysis for Duan's criterion [40] in its actual form or even a general witness that is linear in covariance matrices [42] is a simple generalisation.

Approximate formula for the EPR-variance
Under the assumption that κ ω m γ m motivated in Sec.III C, we made the approximations of Eqs.(11).with ω E/L := ∓(ω m + iΓ) (minus signs associated to index E) and is the normalization factor.The integrands in Eq. (B8) are thus products of the mode functions, Eq. (B13), with the mechanical susceptibility, Eq. (11a), or its modulus square, Eq. (11b).Terms like (ω − ω E ) −1 (ω − ω − ) −1 and (ω − ω L ) −1 (ω − ω + ) −1 are resonant and the next largest terms scale relatively like O Γ n th (C q + 1)/ω m .As a further (purely technical) simplification we keep only the resonant terms and require that Γ ω m / n th (C + 1); this is the origin of Eq. (25).The thus approximated four integrals of Eq. (B8) can be evaluated with the residue theorem, we find [keeping the line structure of Eq. (B8)] (B15) This concludes the derivation of the approximate formula of the EPR-variance, for ideal detection (η = 1), in the regime κ ω m γ m , Γ ω m / n th (C + 1) and ∆ = 0.It is Eq. ( 24) evaluated at η = 1 (called Finally, we model inefficient detection (passive losses) by a beam-splitter of transitivity η right before an ideal (efficiency 1) detector.In the detected channel the field operator is in which c out is the cavity output operator as in Eq. (B5), and c shot is the operator of shot noise that entered the free port of the beam-splitter.The correlators appearing in ∆ EPR , Eq. (B7), factorize because c out and c shot are uncorrelated.Using that the intra-pulse auto-correlators of shot noise equals 2, we find that the EPR-variance accounting for detection inefficiency, Eq. ( 24), is given by We have denoted here with ∆ EPR | η=1 the expression for the EPR-variance for unit detection efficiency derived in this appendix, Eq. (B15).

Figure 3 .
Figure 3. Summary of the scheme we propose: an optomechanical cavity is driven by a continuous, constant and resonant laser (xin, pin).The output field (xout, pout) is continuously monitored by a detection scheme D (e.g., homodyning).ID(t)is the corresponding measurement record (e.g., photocurrent).Its overlap with the appropriate mode functions fE/L, according to Eqs. (17) and (2), yields the quadratures of the early (E) and late (L) temporal modes.

Figure 4 .
Figure 4. Approximate minimal EPR-variance from Eq. (27) (red solid line), exact EPR-variance (green circles) and optimal entanglement witness (blue dots) versus optomechanical coupling g (lower x-axis) and quantum cooperativity Cq (upper x-axis).Threshold for separable states is at 2 (solid black line), the yellow shaded area represents the region of entanglement for the EPR-variance and the optimal witness.This applies to all other figures.

Figure 7
Figure 7. Exact EPR variance (circles) and optimal entanglement witness (dots) versus red detuning of the driving field from cavity resonance.The optomechanical coupling is kept constant corresponding to a quantum cooperativity of Cq = 0.1, 1, 10 for red, green, and blue symbols, respectively.The dashed lines are the approximate minimal EPR-variance formula, Eq. (27), valid only at ∆ = 0 (displayed for comparision).