Enhanced Gravitational Entanglement via Modulated Optomechanics

The role of entanglement in determining the non-classicality of a given interaction has gained significant traction over the last few years. In particular, as the basis for new experimental proposals to test the quantum nature of the gravitational field. Here we show that the rate of gravity mediated entanglement between two otherwise isolated optomechanical systems can be significantly increased by modulating the optomechanical coupling. This is most pronounced for low mass, high frequency systems – convenient for reaching the quantum regime – and can lead to improvements of several orders of magnitude, as well as a broadening of the measurement window. Nevertheless, significant obstacles still remain. In particular, we find that modulations increase decoherence effects at the same rate as the entanglement improvements. This adds to the growing evidence that the constraint on noise (acting on the position d.o.f) depends only on the particle mass, separation, and temperature of the environment and cannot be improved by novel quantum control. Finally, we highlight the close connection between the observation of quantum correlations and the limits of measurement precision derived via the Cram´er-Rao Bound. An immediate consequence is that probing superpositions of the gravitational field places similar demands on detector sensitivity as entanglement verification.


I. INTRODUCTION
One of the central features in most attempts at combining gravity and quantum theory is the assumption that the gravitational field should be quantised.However, in contrast to the other known forces, there has been no experimental evidence to motivate this approach.Instead the reliance is based largely on a combination of aesthetic appeal along with a variety of well known inconsistency arguments for coupling classical and quantum fields [1][2][3].As such, the question of experimental verification, and indeed whether this is in principle even possible, is of significant importance.
In recent years one common approach is to focus on identifying ways that phenomenological features, such as a minimal length or modified dispersion relation, may be accessible in the low energy regime.Surprisingly, this has led to a number of proposals for so-called table-top experiments, in particular within the framework of atomic, molecular and optical (AMO) physics [4][5][6][7][8][9], (see also [10] for a recent review).In these approaches, one attempts to effectively reach the quantum gravity regime through the use of composite systems, provided that large mass quantum states can be produced.
An alternative strategy is to bypass theory specific predictions and instead tackle the question of nonclassicality directly.Recently it has been argued that witnessing the build up of quantum correlations through only a gravitational interaction is sufficient to rule out a purely classical description of gravity [11][12][13].This observation rests on the fact that Local Operations and Classical Communication (LOCC) cannot increase entanglement between two systems [14].From this perspective, entanglement * alexander.plato@uni-rostock.deverification can be viewed as a higher level test of the underlying theory.
One can expect that the detection of quantum features should depend heavily on the sensitivity of the test system to gravity.This requires not only that the field is measurable, but that the precision should be high enough to detect quantum fluctuations of the source.It is well known that by employing resonant driving, measurement sensitivity can be significantly improved.A promising setup is then provided by two disconnected optomechanical systems placed in close proximity.Here, radiation pressure induces oscillations in the mechanical elements, which if sufficiently isolated will only be coupled through their gravitational interaction.Individually, such systems serve as promising candidates for gravity sensors [15][16][17][18], and the fundamental limits on their precision in the nonlinear regime have previously been explored [19,20].An important observation is that a much greater sensitivity can be achieved if, in addition, the optomechanical coupling is modulated at resonance with both the mechanical and driving frequency [21].
This motivates an extension to an optomechanical experiment first suggested by two of us to entangle the motion of two movable mirrors through their mutual gravitational field [13].We show that if gravity can be modeled by a quantised Newtonian interaction then modulations of the optomechanical coupling close to the mechanical resonance can significantly enhance the rate of entanglement generation.Such modulations have been experimentally demonstrated, for example, in nanomechanical setups [22] or in levitated optomechanical systems by employing hybrid-Paul traps [23][24][25].In particular, we find that entanglement can grow asymptotically with the cube of the interaction time, and proportional to the variance of the photon number operator of the initial cavity state (corresponding to roughly the position variance of the mirrors).
In fact, the link between measurement precision and the generation of quantum correlations can be made more concrete and we will show that the integration time needed for an optomechanical sensor to achieve a measurement precision to better than the quantum variance of the corresponding operator (in our case the mirror position or cavity photon number) is closely related to the entanglement time.From this perspective, the development of (decoherence free) quantum systems as sensors is as important as creating large nonclassical states.In particular, this implies that proposals to test the superposition principle for gravitational fields [7,8] are operationally just as difficult as entanglement verification.It also emphasises that one does not always need to resort to strict gravitational cat states in order to test for nonclassical features of gravity.
In this paper we analyse the feasibility of performing such an experiment in a regime where entanglement would be expected.We start by introducing a model Hamiltonian which describes the low energy gravitational interaction of two optomechanical systems in the nonlinear regime.Solving the dynamics for arbitrary time dependent parameters leads to a coupling between the cavity fields proportional only to the respective photon number operators.This allows for a simple estimate of the entanglement rate, provided measurements can made when the combined field-field state is approximately pure.
We then quantify the effects of extraneous couplings through the use of a suitable entanglement witness.This reveals a trade off between the width of the viable measurement window and the initial cooling (previously noted for the interference visibility [26]), but also exposes additional demands on frequency matching of the two mechanical elements.For high photon numbers, these constraints can be significant, though we show that at least for the first, the temperature dependence can be reduced by using purely oscillatory modulations.
On the other hand, decoherence requirements are severe and cannot be mitigated through local control of the dynamics.We recover known noise conditions on the mechanical relaxation rate in the limit of perfect measurement sensitivity [27,28], as well as the more stringent bounds needed to keep the witness within some fraction of its optimal value.In both cases these depend only on the mass, equilibrium separation and the temperature of the environment and suggest improvements of many orders of magnitude are needed before such an experiment can be performed.
For lab based experiments, where the typical energy scale is small, a common approach is to start by quantising the classical Newtonian potential [13,27,29,30], , where G is the gravitational constant and m i are the masses of two particles separated by distance r.
Intuitively this can be done by introducing the quantised particle positions x1 and x2 as small perturbations around an equilibrium distance d, such that r → r = d + x2 − x1 .However a similar result can also be obtained from the gravitational interaction of two scalar fields in perturbative quantum gravity when restricted to the two-particle sector (in the non-relativistic limit) [31] [29].In both cases this leads to a Hamiltonian term, Provided the expected displacements are small, we can expand to second order in x2 − x1 , The last term in particular leads to a quantised gravitational interaction Ĥint = 2Gm 1 m 2 d 3 x1 x2 , and enables the generation of entanglement between the two particles.In general this interaction will not commute with the local dynamics, which in turn leads to non-trivial evolution on the relevant timescale for entanglement generation.The simplest example is given by two free particles with Ĥ0,i = p2 i /2m i , where the additional second order terms in (2) give rise to local Hamiltonians for (shifted) inverted harmonic oscillators, Ĥ = p2 where γ = Gm 1 m 2 /d 3 .
Here the behaviour of entanglement generation can be markedly different if the sign of the potential is reversed [32], and in particular tends to be lower (or indeed can return to zero) when the particles are trapped in harmonic potentials.In practice, however, it can be beneficial to maintain localisation through the introduction of an additional potential (which can be either time dependent [11] or independent [13]).One can then attempt to offset the reduced entanglement rate by transferring the correlations to a set of ancilla states in such a way that they increase after every cycle.The advantage of this approach is that entanglement can often be more easily verified between the ancilla degrees of freedom, either because they are chosen to have a reduced dimension or because they are easier to access in an experiment.This type of setup forms the basis for the majority of proposals to test gravitationally generated entanglement.

III. MODEL
A practical realisation of the above can be provided by optomechanical systems, in which the entanglement is transferred to optical degrees of freedom via a dispersive coupling of the form g(t) N x (where N is the photon number operator).We will now show that for two coupled cavities, the gravity mediated photon-photon coupling takes on a surprisingly simple form: Û12 (t) = e iD(t) N1 N2 .The reduced complexity of the interaction simplifies the analysis of the entanglement rate, and motivates an optimal choice for the initial cavity state.In practice however, this is not easily achieved.Instead we argue that by tailoring the time dependence of the optomechanical coupling g(t), one can significantly improve the performance of less optimal (but easier to produce) cavity states.The moving elements interact via their gravitational fields.To second order this leads to an x1 x2 coupling and allows the generation entanglement between the mechanical modes.This is exchanged back and forth with the cavity fields, which themselves become entangled around multiples of the mechanical period 1/ω m .
We begin by assuming two optomechanical cavities coupled only through the gravitational interaction of their respective mechanical elements.This can be realised, for example, in a Whispering Mode Gallery system where the mechanical elements are made up of fuse-silica microspheres.An alternate setup is shown in Fig 1, using Fabry-Pérot cavities, but in principle the analysis can be readily applied to a range of systems.In general, the precise form of the interaction will depend on the geometry of the masses, m 1 and m 2 , however for simplicity we assume these to be by point sources (2).If each cavity contains only a single mode, then the total Hamiltonian describing this system is given by, where Nj = â † j âj are the number operators for the respective cavity fields (with frequency ω 0,j ), ω2 j = ω 2 m,j − 2γ m j are the mechanical frequencies shifted by the inverted harmonic potential coming from (2), and S j (t) are arbitrary linear displacements.The latter contain the first order gravitational contribution in (2), but could also include effects from additional external potentials.Note that the relative minus sign between cavities accounts for the opposite alignment.The time evolution corresponding to (4) can most easily be solved by first decoupling the mechanical degrees of freedom.This can be achieved via the symplectic transformations [32], where c ≡ cos(a) and s ≡ sin(a), provided that the parameters a and b satisfy, It is convenient to rewrite the transformed mechanical degrees of freedom in terms of the mode operators, b± and b † ± via x± = x 0,± ( b † ± + b± ) (and similarly for p± ), where x 0,+ = 2m 1 ω + , x 0,− = 2m 2 ω − , and, We note for later, that when ω m,1 = ω m,2 ≡ ω m the frequencies simplify to . So far the parameter dependence in each cavity is arbitrary, however we now impose the restriction g j (t) = g(t)ḡ j .This allows us to define time independent combinations of the number operators associated to each plus/minus mode.The Hamiltonian (4) can then be written as Ĥ = 2 j=1 ω 0,j Nj + Ĥ+ + Ĥ− where, (for index µ = +, −), with the transformed number operators, and k ± (t) = x 0,± g(t)/( √ 2 ω ± ).Similarly, the (potentially) time dependent displacement terms are defined by D 1,+ (t) = In general these terms do not contribute to the entanglement, though they can play a role in determining the allowable parameter regime (see section VII).
We now observe that, for j = 1, 2; µ, µ = +, −; µ = µ , and so the the evolution operator corresponding to (4) can be split into, where each Ûµ (t) corresponds to the individual evolution of the Hamiltonians (8).A general ansatz can be found by considering the closed Lie-Algebra generated by the terms in (8) [33], where Dµ (α) = exp(α b † µ − α * bµ ) is the displacement operator acting on the transformed mechanical modes, and the time dependence of the coefficients have been suppressed for notational convenience.Transforming back to the physical cavity field modes, we find the total evolution operator takes the form, 2 ) e iD N1 N2 where, When ( 8) are time-independent, then the solution ( 12) is well known [19,20,34,35].More recently, the extension to arbitrary time dependence has be been obtained in [33].Here the solutions have the following integral representation (note the additional ± notation should not be confused with the transformed modes, µ = +, −, appearing above), where (see [33] for details), For our purposes, the most the relevant terms are contained in the field-field interaction Û12 (t) = e iD(t) N1 N2 , as well as the coefficients K Nµ , which determine the coupling of the fields to the mechanics.The latter leads to a drop in purity of the combined optical state which in turn reduces the observable entanglement.While these terms in general will remain finite at all times t > 0, for certain choices at least one µ-mode can be made to decouple.For example, in the case of a single time-independent optomechanical system, a well known feature is that the field and mechanics disentangle at multiples of the mechanical period, t = 2π/ω m .In our model, this corresponds to, however, a non-zero gravitational interaction means, ω + = ω − , and so both cavity fields cannot simultaneously decouple [36].
In principle this means that the entanglement generated between the fields will always have some dependence on the initial states of the mechanical elements.We will asses the effect of this contribution in section V, however a simple estimate of the entanglement rate can be made by temporarily ignoring the field-mechanics interaction.The evolution operator (13) then has a reduced part, Ûred = Û1 Û2 Û12 , where U 12 is defined as above, acting only on the physical cavity field modes.For interactions of this form, we can show that the linear entropy, S L = 1 − tr[ρ 2  1 ], for an initially separable pure state is given by, (see appendix A for details and the full expression).A characteristic timescale, τ e , can then be found by solving, Thus, the entanglement rate, at least when S L is small, is enhanced by states with large initial variances.For fixed photon numbers, N p ≡ Nj , the optimal choice is a superposition of Fock states, however generating these for large N p is difficult in optomechanical systems.A much higher photon number variance can be achieved by using coherent states, where (∆ Nj ) .Oscillations for large DN p (visible for the blue curve) are a consequence of periodicities in the linear entropy, which begin to appear around the point it reaches its maximum.All plots we made using (A16).

IV. ESTIMATION OF ENTANGLEMENT TIME: CONSTANT VS. MODULATED COUPLING
For (18) to be valid, one needs to evaluate D(t) at a time when both the K µ (t) coefficients are sufficiently small.If the bare mechanical frequencies differ significantly, this may not be possible (we address this issue in section V B 2), and in practice the choice will largely depend on whether one mode is easier to cool.For reasonably symmetric cavities, a sensible approach is to measure at a multiple of the average of the mode decoupling times.In the example above, this corresponds to t q = qπ( where q ∈ Z + .However, in general we will adopt the label t q to correspond to the qth (chronologically occurring) optimal measurement time for a given strategy.

A. Time Independent Systems
In the case of time-independent optomechanical systems, it is straightforward to calculate the relevant coefficient in (12) necessary for estimating the entanglement rate.This is given by [34,35], In the ideal scenario, the mechanical frequencies should be equal, ω m,j ≡ ω m , as this maximises the viable measurement window.We then note from (7), that the difference in mode decoupling periods is approximately ∆t = 4πγ g /ω m .Therefore, provided the number of mechanical periods q 1/γ g , then D(t) will vary slowly between the respective decoupling times.Inserting the above into (14), we find to lowest order in γ g , where are the original coupling constants for each cavity.From (19) we can then estimate that, for coherent states, the characteristic entangling time is given by, In general, achieving coherence for large masses on this timescale is a significant challenge and so the most natural approach is to choose m 1 = m 2 ≡ m for the largest viable m.Realistic parameter values are given in Table I, then assuming k 1 = k 2 ≡ k 0 , suggests an entanglement timescale on the order of τ e ∼ 10 6 s.Thus we see that even with high numbers of photons, the entanglement time is still long in comparison to typical noise timescales.A key result of the present work is that τ e can be further decreased by modulating the parameters of the system.In particular, we will now show that driving the optomechanical coupling at, or close to, resonance with the frequency ω m can lead to cubic time dependence for D(t) in the long time limit.This leads to a significant speedup in the generation of correlations between the optical modes.

B. Entanglement Enhancement Through Modulated Optomechanical Coupling
We now assume that the light matter couplings are modulated according to g j (t) = g 0 + ¯ cos(Ωt + φ 1 ).When  [15], where a low-noise optomechanical system based with a silica mass in the mg range was implemented.R is calculated assuming spherical mechanical elements.For these values γ g ≈ 2 × 10 −18 .
resonant with one of the mode frequencies, this continually drives the (expected) oscillation amplitude, and position variance, of that mode [21].As a result, entanglement generation is increased through the field-field channel.However, unlike in the constant coupling scenario, the mechanical state will vary after each oscillation and so no times t > 0 exist when the field and mechanics factorise, i.e K µ (t > 0) = 0.The solution is to instead modulate the coupling at specific fractional frequencies of the form Ω n,µ = (1 − 1/n)ω µ , where n ∈ Z.It was shown in [21] that for times t q,µ = 2nqπ/ω µ the cavity and mechanical states do indeed disentangle.Evaluating the F N 2 µ , and F Nµ Bµ,± integrals in (16) we find that for an optimal phase choice, φ 1 = π/2, and for times t q,− < t q < t q,+ , the D coefficient is given by, to first order in γ g = Gm/ω 2 m d 3 , where k 0 = x 0 g 0 /( ω m ) and = x 0 ¯ /( ω m ).Note to this order, the choice between Ω n,± is unimportant here.
While D(t q ) grows linearly in q, the term corresponding to modulated optomechanical coupling grows cubically in n.This suggests the most effective strategy is to measure at the first optimal measurement time (q = 1).It is then clear that for ω m t 1 1, the second term will dominate (provided is not much smaller than k 0 ), that is, when multiple mechanical periods are needed to entangle the systems.For optomechanical cavities ω m is typically larger than 2π × 10 rads −1 , and so if entanglement can be reached within t ∼ 10 −2 s, one would choose to work with unmodulated systems.The reality however, is that even for optimistic parameter choices, k 2 0 γ g N p 1, and so ω m t 1 needs to be large.This leads to our main result: in the long time limit, |D(t 1 )| ≈ 1 4π 2 2 ω 3 m γ g t 3 1 , and so the entanglement time for cavities initialised in coherent states is given by, By way of comparison, using the parameters in Table I and setting = k 0 , modulation decreases the entanglement time significantly to τ e ≈ 2 s, which corresponds to an improvement by 6 orders of magnitude.In general, the modulated entanglement rate is enhanced by a factor of 1 12 ( /k 3 0 ) 2/3 (πN p γ g ) −2/3 compared to the constant coupling case.This means it is particularly advantageous in precisely the experimental regimes which are most realistic: weak gravitational interaction, small light-matter coupling and low photon numbers.
In practice, the restriction to fractional frequencies does not significantly decrease the direct field-field coupling over the resonant case.Similar calculations for, e.g.Ω = Ω ∞,+ = ω m , lead to, where we have taken t q = 2qπ/ω m .For ω m t q 1, the last term in equation ( 25) dominates and so compared to (23) at a given target time (choosing the driving frequency as Ω n=q,± ) we see that the modulus differs by only a factor π 2 /3 ∼ 3. On the other hand, the larger K N µ (2qπ/ω m ) values resulting from the resonant coupling impose stricter cooling requirements for the mechanical modes.

V. ENTANGLEMENT VERIFICATION AND NOISE
The results of the previous sections correspond to the best case entanglement rate between the two cavity fields.In reality, a combination of external decoherence along with other intra-system couplings will mean that the final field-field state will not be pure.In general, quantification (and measurement) of entanglement in mixed states is typically challenging, though there are notable exceptions -as an example, in Appendix D we analyse the logarithmic negativity for initial cavity states of the form |ψ = 1 √ 2 (|0 + |N ), recovering the qualitative features outlined in the following sections.Instead it is often more practical to resort to so-called entanglement witnesses, Ŵ, which verify with certainty the presence of entanglement if W = Tr Ŵρ < 0, but can make no statement otherwise [37].
From a technical perspective, this means one must be careful when using such operators to infer an entanglement rate.In particular, Tr Ŵρ(t 2 ) < Tr Ŵρ(t 1 ) does not guarantee that entanglement has increased in time.On the other hand, achieving W(t) < 0 to begin with may impose a minimum constraint on the allowable noise in the system, and this can have a meaningful time dependent operational interpretation.We will quantify this in the sections below.
In practice, however, one will also need a sufficiently high signal-to-noise ratio [38] in the measurement of W on the state for entanglement to be observable.The particular threshold needed (i.e.W < −thresh.) will depend heavily on the experimental apparatus used.As such, identifying this cut-off a priori is difficult, but we can maximise the chances of verification by performing the measurement when the witness function reaches its minimum value below zero.
In the following section we identify a witness that not only detects entanglement in our system, but also recovers approximately the same timescale (24) in the limit that all field-mechanics interactions vanish (at the chosen measurement time).More generally, this scaling is preserved as long as the noise constraints (which we will derive explicitly) are satisfied to within roughly an order of magnitude.This holds both in the case of extraneous intra-system couplings and for thermal decoherence.

A. Entanglement Witness
In order to identify an appropriate entanglement witness it is useful to anticipate the ideal joint state of the field modes.Intuitively, when the gravitational coupling strength is large, then the mirrors effectively become rigidly connected and so the system is analogous to a single cavity containing two modes [13][39].For certain parameter choices, it is known that these modes can become entangled, forming the two-mode cat state, [35].However, many of the simplest continuous variable entanglement criteria fail to witness entanglement in these types of states.One of the first examples that does was provided by Shchuckin and Vogel [40], and can be repre-sented by the following determinant, This can be measured, in principle, using homodyne correlation measurements [41].
In Appendix C we provide explicit evaluations of the time dependent terms appearing in ( 26) assuming an initial state of the form ρ(0 , where |α is a coherent state and ρ th µ are thermal states of the mechanical modes.For simplicity, we will continue to work almost exclusively in the symmetric cavity approximation (though the extension to arbitrary parameters will turn out to be straightforward).Under these conditions cos(a) = sin(a) = 1 In general, the resulting expectation values are functions of the absolute values, |α 1 | and |α 2 | of the coherent state parameters for the two cavity modes, and not the relative phase.In the symmetric cavity scenario a natural choice is to set these equal, α 1 = α 2 = α, with |α| 2 = N p .Even then, the full witness expression is still difficult to treat analytically, however, numerical results suggest that the restriction to B = 0 at the measurement time provides a somewhat optimal regime.While this term has no effect on the entanglement, similar terms are known to influence the sensitivity of optomechanical sensors under certain POVMs (e.g.homodyne measurements, [21]).For = 0, it is usually sufficient to consider only integer values of k 0 , but more care must be taken in the modulated scenario.In these cases, the witness is given by the simpler expression (C25), which holds both when κ µ describes the effects of the traced out thermal mechanical modes and, as we will see later, decoherence from an environment.It is therefore convenient to write κ µ = κ th µ + κ dec µ , where in the present section we will focus on the first term, given by κ th µ = |K Nµ | 2 (n µ + 1/2), where nµ = 1/(e ω s /k B T µ − 1) is the occupation number for temperature T µ (see (C8)).
In general, the optimal measurement time will lie some- where between what would be the individual mode decouplings.For equal temperatures, this turns out to be just the mean, t q = qnπ(1/ω c + 1/ω s ), where to preserve symmetry we will from now on choose the modulation frequency as To lowest order in γ g , and for large n, both κ th µ (t q ) are equal to (C13), Neglecting the decoherence contribution, we can then set κ µ = κ in the witness (28), which further simplifies to where y = 2N p z 2 + κ(t q ) and z = sin(D(t q )/2).
We can now apply the requirement that W 1 (t q ) < 0 to give an upper bound on the thermal noise term.This leads to, If an experiment is sensitive to small amounts of entanglement, such that very early measurements are viable (D(t q ) 1/N p ) then we can expand the right hand side to first order in z, This gives the lowest acceptable noise for which entanglement can be verified.We find a similar condition for the logarithmic negativity when the initial cavity fields are in superpositions of Fock states (see (D7)), suggesting the behaviour of the witness is reasonably faithful for small amounts of entanglement.In particular, for κ = 0, entanglement can be immediately detected for any non-zero D(t q ).On its own, however, this is not especially useful, at least in the case of purely thermal noise.This is because κ th (t q ) and D(t q ) have different time dependencies, and indeed the former tends to grow faster than the fieldfield coupling.As a result, (32) only provides an upper bound on when entanglement can be detected.To see this explicitly, it is convenient to separately consider the constant and modulated-dominant regimes.For later, we will also find it useful to regard κ in general as a function of z ≈ D(t q )/2, then using ( 29) and (23) we have to lowest order in γ g .Applying the inequality (32), and inverting for t q we find that in both cases (up to a small numerical difference), t q 12/(ω m γ g (2n + 1)).Thus, thermal noise limits how long an experiment can be performed.In fact, a finite measurement time is a generic feature of our chosen witness -even for κ = 0, the next order expansion of (31) leads to |D max (t q )| ∼ 2|z| 2 1+2N p , despite the fact that the maximum entanglement has not yet been reached (see figure 2).
A more useful indicator for the ideal measurement time is to instead look for the point of maximum witness violation.This can seen in figure 3, where W 1 (t q ) reaches a minimum of − N p 4 in the limit κ th → 0. We can then use this as a reference point to determine the maximum thermal noise for which the witness only drops to some fraction of this value.

B. Optimal Measurement Time
To proceed, it will be convenient to take κ = κ(z), then the minimum of ( 30) is achieved for z min satisfying, 4z 2  min + e 2y min = 1 + If κ varies slowly with respect to z (i.e.κ (z) 4N p z) we can ignore the noise gradient term in the denominator.Substituting for y min we then have, Now, for realistic values of γ g , z min must be small on accessible timescales.This suggests N p 1 (which we will assume from now on), and so expanding the logarithm to first order, we have after a simple rearrangement [42], Thus, in the κ → 0 limit, the witness will take its minimal value when D(t q ) ≈ 1/ √ 2N p .Substituting into (30), we have, as claimed.At high temperatures, the slowly varying condition is typically not satisfied, and so the derivative term in (35) must be included.Following from (33), we can write κ th = Γ th 0 z 2 and κ th = Γz 4/3 , for the constant and modulated optomechanical coupling regimes separately.The minimum conditions in the first case are easily found by noting that (34) is equivalent to (35) with the substitutions N p → N p + Γ th 0 /2 and κ → 0. We then find, which implies the measurement time should satisfy, (where we have again used N p 1).Thus, just as with D max , the optimal timescale shifts earlier at higher temperatures.This is, of course, not an advantage because the added noise leads to an increase in the minimum value of the witness (see figure 3).At the first optimal measurement time, this is given by, In effect, it prevents the experiment from running long enough to build up a highly entangled state.We can characterise an acceptable level of noise by demanding that the magnitude of the witness decreases by no more than some factor, a, of the idealised (κ th = 0) value.This amounts to solving W 1,min = aW 1,min | κ→0 , for which we Therefore, to keep the witness below, for example, half its absolute minimum (for a given set of experimental parameters), the mechanical modes must both be cooled to, As N p will typically be very large, this should not pose a significant obstacle.
For modulated optomechanical couplings, the analysis proceeds along similar lines though it is convenient to first make a suitable expansion of (34) (see appendix C for details).In this case we find (to first order in 1/N p ), 3/2 , and where v = 5 18 It is immediately clear that for all but very large n, the minimal witness time is very close to that for κ = 0. Indeed, substituting into (30), and choosing a = 1/2, eventually leads to the bound, nmax 0.66 Thus, for realistic experimental parameters (where γ g 1) modulated couplings significantly improve tolerance to thermal noise.Indeed, for the values in table I and = k 0 , n 4 × 10 13 , which corresponds to a mode temperature of T ≤ 2 × 10 6 K.

Thermal Dependence of the Measurement Window
The constraints above are valid when the measurements are localised to within a very short window around t q (which we can expect is less than the decoupling time difference, δt ∼ 2γ g t q ).It is important, however, to know the extent of the usable measurement window, both because such processes are not instantaneous and because reducing the signal to noise ratio will rely on obtaining long integration times.The easiest way to estimate this is via the condition (31), along with expressions for D(t q ±δt) and κ th (t q ± δt).
This approach will only be viable provided B(t 1 ± δt) ≈ 0, otherwise the full witness expression must be used.In general, however, the window will be small on the timescale of the mean mechanical period τm = π(1/ω c + 1/ω s ).Expanding in terms of δt = rτ m and assuming γ g < r 1, it can be shown that, (where additional terms only enter at order γ 2 g r 2 ).Similarly we can find, and, The first of these expressions tells us that even at zero temperature, B will be much smaller than κ th , and so we can reasonably expect (28) to hold.We also note that for κ th the variation in r does not depend on γ g , and so thermal noise grows much faster than the coupling D. This means the measurement window will be much smaller than τm (as D(t q ) ∼ O[γ g ]), and we therefore take D to be effectively constant over this timescale.Moreover, when n nmax , we can also ignore the first term in (45) and so around the decoupling time, κ th is independent of D. This immediately allows us to determine both the conditions and size of the maximum measurement window: from (31), the maximum κ th for which the witness is exactly zero occurs around t q satisfying, where for small D, we take the positive solution, and so, with a corresponding κ th max (t q + rτ m ) ≈ 1 8N p .Substituting in (45) when n nmax , we find that the maximal size of the measurement window is [43], At zero temperature, this equates to around ±36ns for the system outlined in table I, which is approaching the bandwidth limits of state of the art detectors.We note that in this case the window size is dominated by the constant component of the optomechanical coupling.If one can instead engineer k 0 = 0, > 0, then the next nonzero contribution to κ th is 4th order in r, and so the maximum δt is instead given by, By comparison, the zero temperature window is increased to ±3.4µs, while making no change to the measurement time.
More generally, we can estimate the measurement window at an arbitrary decoupling time by evaluating (31) using (45) and D(t q ).In the physically relevant regime, it is enough to expand (31) to second order in D(t q ).We then have, where again we have taken N p 1.At the witness minimum, we typically have D(t q ) ∼ 1/ √ 2N p , and so the window at this time is still around 91% of its maximum value.Surprisingly, the estimates ( 49) and ( 52) appear to work well even at temperatures where the analysis breaks down, see figure 3.

Frequency mismatch
The results of the previous sections are reasonably tolerant to most sources of asymmetry.However, the most obvious exception is when there are mismatches in the mechanical frequencies.This is because the separation between mode periods grows and so it becomes much harder to find a time when both κ th µ are small.As the analysis in this case is significantly more complicated, we estimate the dependence only in the unmodulated scenario.Using the general frequencies ω µ (7), we insert the B µ coefficients from ( 20) into ( 14) and expand D(t) to lowest order in the (sum of the) relative frequency shifts in ωµ , i.e.
wlog), we find after some work, where k j = x 0,j g 0 ḡj ω m,j , with g(t) = g 0 , are the original coupling constants for each cavity.In deriving the above FIG.3. The (rescaled) witness as a function of discrete t q times for constant optomechanical coupling (n = 1), at both n = 0 (a) and the predicted nmax given by (41) (b).Here we choose the more optimistic parameters, ω m = 2π × 10 2 Hz and m = 2.4 × 10 −8 kg, with N p = 10 6 , leading to γ g = 5 × 10 −13 .Between any two t q 's, the witness can only verify entanglement within a small measurement window |t − t q | δt (see insets).Numerical results agrees well with the theoretical value (52), even beyond its apparent regime of validity.However, at high temperatures, the t q with maximal measurement window no longer follows (48).Note, there exists a maximum verification time (see (68)), after which the witness is effectively always positive, even though entanglement may still be expected.
it is useful to adopt k + = ( 2 ).We now assume a small frequency mismatch, ω m,1 = (1 + r f )ω m,2 , where r f is a positive constant.Around a measurement time t q = πq(1/ω + + 1/ω − ), we have, where D(t q ) is as given in ( 21) (for ω m,2 = ω m ).Provided r f 1, we again see that D(t) varies slowly around t q .A similar analysis can be performed for κ th µ terms, making use of (17).For r f much larger than the relative frequency shifts, we find where again we are considering only the nµ nmax limit so that κ th µ (t q ) can be neglected in comparison to D(t q ).Note that the expressions above do not depend on the gravitational field, which could have already anticipated from (45).Thus τm can effectively be taken as the mean of the bare mechanical periods, which makes the calculation significantly easier.For the sake of simplicity, we assume k j = k 0 , nµ = n and r ∼ r f then in the long time limit (large q), κ th + ≈ κ th − is dominated by the r f term.Using the analysis of the previous section, we immediately find, Thus we see the trade-off between controlling thermal noise and the matching of mechanical frequencies -even for perfectly localised measurement times, a non-zero ∆ω m can severely limit the upper bound on n.Note, the difficulties highlighted in the previous two sections are not inherent to optomechanical systems alone.Recall that the noise term κ th (t) is directly proportional to |K Nµ (t)| 2 , which itself is related to the total displacement in phase space (as can be seen from ( 13)).Thus the measurement window acts as a proxy for the level of control needed to return the mechanical degrees of freedom to close to their initial state.

C. Thermal Decoherence
The relative insensitivity to the initial mechanical state arises because K Ns (t 1 ) can still be very close to zero even at the centre of mass mode decoupling time (or viceversa).A more significant experimental challenge comes from external (open-system) noise processes, which in general cannot be avoided by choice of the measurement time.These can be broadly split into those acting directly on either the optical or mechanical degrees of freedom [44].Optical losses, for example, can be expected to reduce the visibility of the signal [45], while mechanical decoherence will prevent the initial build up of entanglement between the masses.
The latter in particular is well known for limiting the generation of macroscopic (spatial) superpositions.In a typical analysis, one modifies the usual von Neumann evolution equation by the inclusion of additional superoperators describing the effects of the environment (along with a corresponding renormalisation of the system Hamiltonian).If the temperature of the environment is assumed to be high (with an ohmic spectral density) then it is usually sufficient to consider only the decoherence term, L In this limit the coefficient is time independent, and given explicitly by Υ = 2mk B T γ R where γ R is the relaxation rate and T is the temperature of the environment [46].The extension to multiple mechanical modes is straightforward as long as the intra-system couplings are weak.In this case a good approximation to the dynamics can be found by adding additional L[ρ] terms for each system degree of freedom [47].Neglecting optical noise, our optomechanical model can then be described by the following equation, where Ĥren is the renormalised system Hamiltonian (which for simplicity we will set equal to H, as this typically amounts to a shift in the mechanical frequencies), and γ R = ω m /Q is the relaxation rate for an optomechanical cavity with quality factor Q.
In order to evaluate the witness (26), we only need to solve for the reduced total optical state.For the single cavity (time-independent) setting, a convenient approach has been presented in [48] where the analogue of the master equation ( 57) is first projected out in terms of the operator ρ D (t) = p, q|ρ(t)|m, n , where |m, n were Fock states of two arms of an interferometer.This can be solved exactly when making use of the partial trace over the mechanical modes, leading directly to the reduced field-field density matrix ρ cav p,q,m,n (t) = Tr M [ρ D (t)].The extension to two cavities is largely straightforward once we move to the decoupled mode picture, in particular where we note that for identical cavities, . For completeness we provide full details in appendix B, where an explicit solution is given by (B20).This can be compared directly to the noiseless cavity reduced state (D2), where we find that the modifications can be completely described by formally replacing the initial thermal occupation numbers, nµ , at time t, by where As these terms always appear multiplied by an additional |K Nµ | 2 factor, (such that the thermal decay coefficient κ µ → κ th µ + κ dec µ ) any potential divergence in the second term is of no concern.Hence, while K Nµ (t) can in principle be zero, κ dec µ (t) is not in general.We also stress that nµ (t) should not be interpreted as the current temperatures of the mechanical modes, but rather only as an effective model for the decoherence induced loss of purity of the cavity state (which occurs even for perfect decoupling).
At zeroth order in γ g , and for n nmax we find that for identical cavities, κ dec µ (t q ) = κ dec (t q ) κ th (t q ), where Thus, from (32) one immediately sees that entanglement can only be witnessed when, which is independent of both the local dynamics and N p .A similar behaviour is found for initial cavity states of the form |ψ = 1 √ 2 (|0 + |N ), where |0(N ) are Fock states.In this case entanglement can be calculated directly using the logarithmic negativity (see Appendix D), which shows this result is not unique to our chosen witness.In fact an identical expression has also been found in the linearised regime for Gaussian initial states [28] (in the constant coupling scenario), where it was conjectured that the bound is a generic consequence of the master equation (57).Thus the 'size' of the superposition created does not matter as both the entanglement and decoherence rates are enhanced in the same way.Note, we have not considered any decoherence on the optical modes.That is, "shifting" the entanglement to an apparently protected degree of freedom does not help, as the mechanism used to transfer the correlations also provides a channel for the decoherence to act.
In practice, however, other noise timescales will mean it is still highly desirable to reduce the entanglement/verification time.Thus a large photon number, or more generally photon number variance, is advantageous.On the other hand, the cavities themselves, along with their relative positioning, limits the extent to which the mechanical elements can be driven.These restriction will be discussed in section VII.
The analysis for the optimal witness time is straightforward, and we leave the details for Appendix C 1. Writing κ dec (z(t q )) = Γ dec z(t q ), where from (60) Γ dec = 2Υx 2 0,m 2 ω m γ g , the minimum is found to occur at, (62) Similarly, the condition for the corresponding witness value to be at least a times the theoretical minimum (κ = 0) is given by, Choosing a = 1/2 then requires Γ dec 0.219, and so γ R k B T 0.11 Gm/d 3 .As before, increasing the noise shifts the minimum earlier in time until above Γ dec = 1 when no entanglement can be verified.
For the cavities listed in Table I satisfying the bound (61) requires a quality factor Q ∼ 10 23 .This is some 15 orders of magnitude higher than what has currently been achieved (see for example [15]), and poses a severe obstacle to witnessing gravitational entanglement.

Measurement Window
Given the close connection between κ dec (t q ) and D(t q ), one may also suspect that κ dec varies slowly around the measurement time t q .Indeed, we find where the dependence on the constant (i.e.k 0 ) component enters at order γ 2 g r 2 .This means the measurement window is dominated by the initial temperature of the state, even if n is very small.We can then substitute into (31), where the second term arises from thermal contribution ( 45), and we have assumed κ th (t q ) is small (i.e.n nmax ).Expanding to second order leads to, which generalises (52).

VI. THE CASIMIR EFFECT AND STRAY CHARGES
In order to maximise the gravitational interaction, the mean separation between mechanical elements should be kept as small as possible.This gives the best chance of satisfying (61).On the other hand, at small distance scales the Casimir force [49] becomes significant and can readily dominate over gravity.A direct comparison can be made via the dimensionless coupling parameter γ C := 1 2 ∂ d F c /mω 2 m , which should be small compared to γ g .In general, the Casimir interaction will depend strongly on both the geometry and material properties of the mechanical components.However when the masses are assumed to be spheres, a straightforward estimate can be made in the zero temperature regime, provided the separation is much larger than the radii, d R. Using equation ( 21) of [50], one finds that, The requirement that γ g γ C immediately implies m 161 π R d 3 m p , where m p = c/G ≈ 2.176 × 10 −8 kg is the Planck mass.This suggests that the masses in general should be (relatively) large, otherwise the separation must increase to compensate.Indeed, for the parameters in Table I, this leads to γ C ∼ 1 × 10 −17 and γ g ∼ 2 × 10 −18 , and so Casimir mediated entanglement would dominate over gravitational effects.On the other hand, a large mass will also be harder to cool which can lead to shortening of the available measurement window, c.f. (66).
To get around this problem, we can adopt an approach frequently used in small distance tests of gravity.This involves introducing a suitable shielding material, for example gold [51], between a driven source and detector.If the material is rigid enough, then the only oscillating signal is that from the gravitational interaction, which can be identified over the stronger (constant) Casimir/electrostatic background.A similar screening effect has been proposed for entanglement generation [13,52], where the transmission of quantum correlations are effectively suppressed because any induced motion in the shield, including its position variance, will be small.In this sense, we can consider the shield as restricting non-gravitational forces to only a classical channel.
It is worth pointing out that the witness ( 26) itself also provides some limited protection against false positives.We note that our previous results can be immediately generalised to γ g → γ tot , where now accounts for all x1 x2 couplings in the system.For example, electrostatic interactions can be included via the Coulomb force, F Q , leading directly to the parameter , where q 1 and q 2 are charges on the mechanical oscillators, and ε 0 is the vacuum permittivity.Thus, in the absence of significant shielding, we can write γ tot = γ g + γ C + γ Q .Now as D(t q ) ∝ γ tot , if one of the additional interaction parameters is comparable to γ g , then D(t q ) will double, leading to a higher than expected entanglement rate.However, as can be seen from figure 3, beyond a certain time the witness is positive.We can estimate this from (31) along with the relevant κ(t q ).For example, in the case of constant optomechanical coupling (33) (and neglecting thermal decoherence), a second order expansion in z leads to the condition, This tells us that above |D(t q )| ∼ 1/(N p + Γ th 0 /2) the witness is positive, and so from (38), if γ tot is larger than expected by ∼ √ 2 an experiment aiming to measure at the optimal time will fail to detect entanglement, even if it is present.For the parameters in table I this would already happen for only 0.006 stray electrons on each mass (ignoring the Casimir interaction).Of course, as z(t) = sin (D(t)/2), there are further periodic solutions, for example when D(t q ) is shifted to D(t q ) = 2π ± D(t q ).However, in this case γ tot would need to be larger by a factor ∼ N p .It then also becomes much harder to satisfy κ < z (from (31)), as κ th (t q ) grows quadratically in γ tot .Thus, in principle, it should be possible to engineer an experiment such that if entanglement is witnessed, we can be reasonably confident it did not emerge from an undesired interaction.

VII. LIMITATIONS DUE TO THE MECHANICAL MOTION
In equation ( 4), we started from a Hamiltonian description of the light-matter coupling based on a perturbative treatment of the motion of the mechanics [53][54][55].This means our results can only be seen as reliable if the motion of the mechanical elements is small.For example, in Fabry-Pérot cavities the displacement of the end mirror must be much smaller than the cavity length, while in levitated systems the relevant lengthscale is set by the wavelength of the light modes [56].Similarly, the gravitational interaction between the mirrors was described by a second order expansion of the Newtonian potential, which requires their displacements to be much smaller than their mean separation.In principle, higher order treatments can be possible (see, for example [57]), however the dynamics will always be limited by the physical dimensions of the experiment.
As the gravitational interaction is weak, by far the strongest effects on the mechanics comes from the radiation pressure of the cavity fields.This shifts not only the equilibrium position, but also increases the oscillation amplitudes of the mechanical elements.The average radiation pressure force, however, can be cancelled using an appropriately tailored external potential (which enters through the C µ coefficients, ( 16)).Thus the mean position can be set to zero (see (E1)), and all that remains are the quantum fluctuations of the mechanical motion.The corresponding variances for the center of mass mode and the stretch mode are (see appendix E) where is the co-rotated position operator of the mode k and for a thermal state (∆x ω µ ) 2 = (2n µ + 1)x 2 0,µ .Imposing the condition ∆x c , ∆x s d (and similarly for the cavity length L) implies a limit on the thermal occupation and the photon number variance.The former has to be minimized in any case to detect entanglement generation, however it is generally advantageous to maximise the latter and so its restriction needs to be accounted for when choosing system parameters.

VIII. COMPARISON TO SEMI-CLASSICAL ANALYSIS
For ancilla based entanglement tests, an alternative approach for analysing the entanglement rate is to work in a semi-classical approximation.This has the advantage that the full, coupled, dynamics need not be solved (in fact, an explicit Hamiltonian is not even required).In the case of optomechanical systems, this turns out to be remarkably accurate when only second order perturbations are included.
We start by constructing a model for two massive particles, where each exists in superpositions of approximately classical trajectories x 1,n (t) and x 2,n (t).These paths are labelled by sets of orthonormal quantum states, |u n and |v n , respectively, encoded on two additional (external or internal) degrees of freedom.We can write the total state of each system as, where |Φ i,n are the mechanical states corresponding to the trajectories x i,n (t) = xi |Φ i,n (t) .The gravitational interaction between the two particles can be modelled by assuming the potential energy difference between each localised superposition gives rise to an accumulating phase, and so after some time, the total state of the combined system is, (74) Note that entanglement only enters at the second order expansion (and above) of (73).We now make the assumption that after a period t q the mechanical and label states disentangle, such that the reduced state of the total ancilla system is pure, For low dimensional systems, the entanglement can be calculated directly from the above.However, the phases φ nm can also be viewed as the eigenvalues of a suitable operator with eigenstates |u n |v m .The former is dependent on the difference between classical trajectories, which suggests a natural construction is to first define two local operators, (76) If the paths can be written as functions of the label eigenvalues u n , and v m , then the desired operator can be found from ∆x(t, Ô1 , Ô2 ), where ∆x(t, u n , v m ) := x 2,m (t) − x 1,n (t).Thus the time evolution of label states (at the decoupling time) is of the form, Û = Û1 Û2 e i Gm 2 t 0 ∆x(t, Ô1 , Ô2 ) By expanding ∆x(t, Ô1 , Ô2 ) −1 to lowest order such that the time evolution operator can be expressed as we can make direct use of ( 18) and the results of Appendix A.
Example 1: Interferometry.One of the simplest examples has been presented in [11], where each particle is instantaneously driven into a cat-state superposition with separation ι, and then returned to its original state after time t.We can define label states |u ± = |v ± ≡ |± ι 2 , then the interaction operator is given by, where Ô1 and Ô2 are constructed from (76).For a two level system, the linear entropy is a function of only ∆ Ô1 ∆ Ô2 = ι 2 4 and so from (94) the condition for maximal entanglement is [11], A characteristic time scale for more general initial states (at least to reach small levels of entanglement) can easily be found from the analogue of ( 19) or, to higher order (A1).
Example 2: Optomechanical systems.In this case the difference between classical trajectories is given by ∆x(t, u n , v m ) = d + √ 2x s (t, u n , v m ), which for thermal mechanical states can be readily found from (see Appendix E), As this depends on the photon number of the ancilla (cavity) state, we set Ôi = Ni .A second order expansion of ∆x(t, Ô1 , Ô2 ) −1 then leads to the time evolution operator where we can identify For both modulated and unmodulated optomechanical couplings we find that D(t q ) evaluates to ( 21) and ( 23) to first order in γ g .The entanglement rate estimates of section IV then follow immediately.Note however, that at higher order in γ g , deviations between the fully quantum and semi-classical treatments begin to appear.

IX. QUANTUM METROLOGY
We can gain further intuition for the observed initial state dependence by making use of techniques from quantum metrology.In this framework, we consider the two cavities as a sensor-source pair.Then a prerequisite for any experiment is that its sensitivity to gravity is high enough to resolve quantum features in the source.In our case, fluctuations of the photon number.
Formalising this perspective is complicated by the fact that there is no parameter corresponding superposition size in the evolution operator (13), even within the semiclassical model discussed in section VIII.This prevents a direct application of the usual methods from metrology.One solution is to first characterise the measurement precision for the case that |ψ source (0) is an eigenstate of one of the quantum degrees of freedom, i.e.Nsource .Then we can effectively replace the operator by its eigenvalue, λ N , in Û and treat this as a parameter for estimation.The quantum Fisher information (QFI) of the evolved state then sets the ultimate lower limit on the variance of an unbiased estimator λN via the Cramér-Rao bound, where N is the number of measurements.We then argue that resolving quantum fluctuations in the source photon number is only possible when ∆ Nsource ≥ ∆ λN .
The QFI is most easily evaluated when the cavity field states are pure.While this will strictly never be true, it is comparable to the approximations used in section IV where we assumed K Nµ (t q ) ≈ 0. In this case it takes the simple form, where |ψ λ N = Û |ψ(0) is the field state of the detector at t q , and |ψ λ N is its derivative with respect to the parameter λ N .Denoting N2 = Ndetector and N1 → λ N in ( 13), a straightforward calculation leads to, (neglecting K Nµ contributions).This is a direct consequence of the approximation, Û12 (t q ) = Û1 Û2 e iD(t q ) N1 N2 .It then follows from (83), and the condition ∆ Nsource ≥ ∆ λN , that Thus we see that metrology predicts a timescale that is broadly in line with that needed for entanglement generation (19).This should come as no surprise as the semiclassical model above tells us that entanglement is only possible when (at least some of) the phases φ nm are different.On the other hand, to distinguish between them, we should be able to resolve differences in the corresponding trajectories, which in turn depend on fluctuations of the photon number.
To make this link more apparent, we can instead derive bounds explicitly for the mirror position.The result turns out to be essentially the same whether we use expectation values or eigenvalues to identify a suitable estimator.Note, while our focus has been on entanglement generation, we could equally interpret the sensitivity requirements as those needed to test whether gravity respects the superposition principle [58].The latter can be traced back to proposals by Feynman [59], but has also received renewed attention in recent years [7,8,30].
To begin, we consider measuring the classical position x source of a particle (of mass m) using a single optomechanical system.In the rotating frame, the detector Hamiltonian takes the form (8), (µ → detector) where by analogy with (3), the displacement due to the gravitational interaction is given by It has been shown in [21] and [60] that the QFI of the local (pure) field state, with respect to the parameter x is given by where (for notational convenience we abbreviate source and detector labels to 1 and 2 respectively).The F -coefficients can be determined by substituting (87) into formally identical equations to (16).Note, in our case the position (as an estimation parameter) has a time dependence, t, which is independent of the actual measurement time, t q .Indeed, if x source (t) ≡ x1 (t) , then for a thermal state (see appendix E1) which can also be zero.Thus, it is convenient to adopt the new parameter N 1 ≡ N1 and write W x 1 = (∂N 1 /∂x 1 )W N 1 .The calculation can be simplified further if we ignore any displacement of the source due to the detector (i.e.setting D 1 (t)) = 0), which means C 1 (t) = 0. Then evaluating W N 1 when both source and detector systems are identical, we find to first order in γ g , which is not a function of the time t.Thus the Cramér-Rao bound for an unbiased estimate of x source (t) using measurements at decoupling time t q = 2πnq/ω m is given by, We now argue that superpositions of the gravitational field can only be resolved when the detector precision, (92), is smaller than the variance of the position operator xsource , (93) where xω (t) = xsource (0) cos ω m t + 1 mω m psource (0) sin ω m t.If we are interested in the entanglement between cavity fields, then the important correlations are those arising between the respective number operators.We can therefore ignore the second term above, and the condition ∆x source (t) ≥ ∆x source (t) immediately recovers (86).

X. CONCLUSIONS
The application of the LOCC framework to tests of quantum gravity has gained significant traction over the last few years.In particular, because direct access to the state of the gravitational field -for example, by probing the existence of gravitons -is almost certainly unachievable.Instead, it is hoped that correlations arising from the collective interactions of macroscopic-scale quantised masses may provide indirect evidence for the nature of gravity.
Such tests are exceptionally difficult, requiring highly non-classical states or very long integration times.This poses a challenge for traditional optomechanical implementations, where generating large cat-states requires superpositions of very different photon numbers.Instead one must typically resort to a combination of low mechanical frequencies and high coupling strengths.
In contrast we have shown that the entanglement rate can be significantly increased if the optomechanical coupling can instead be modulated close to resonance [61].For k 0 = , this advantage scales as roughly -i.e. the benefit is greatest precisely in the most accessible parameter regime.This can be traced back to the position response of the differential mode, which to a good approximation governs the coefficient D(t), and not the initial cavity state.The latter can of course be leveraged for further improvements.For example, (19) can also be readily applied to squeezed coherent states [62], where for squeezing amplitude r and coherent state parameter α the photon number and its variance are given by N p = |α| 2 e 2r + sinh 2 (r) and (∆ Nj ) 2 = |α| 2 e 4r + 1 2 sinh 2 (2r), respectively (for appropriately chosen phases) [21].For squeezing of around 10 dB [63], corresponding to r ∼ 1.73, the requirement on D(τ e ) would be decreased by an additional factor of 30.
In the absence of noise, modulations would therefore suggest that entanglement timescales comparable to proposals based on mechanical cat-states, (e.g.[11]) may be realistic.However, even with these improvements, significant obstacles remain.First, the measurement window for witnessing entanglement is extremely sensitive to the initial temperature of the mechanical modes.Similar behaviour is well known for the single photon interference visibility width in a Michelson style optomechanical interferometer [26,64], which scales roughly as 1/n.In our case, the additional photon number dependence in ( 49) is particularly problematic if entanglement is to be achieved on short timescales.To some extent, this can be alleviated through purely oscillatory couplings (i.e.k 0 = 0), which reduces the dependence to ∼ 1/(N p n) 1/4 .However, in practice this still means that the system must be cooled far below the limits set by (43), even for state of the art detector bandwidths.
Secondly, in a two cavity system, the relative dynamics must be matched to extremely high precision.This guarantees that both mirrors simultaneously disentangle from the optical modes (at least approximately), maximising the measurable entanglement.Provided both cavities are initialised at the same time, the main determiner is whether the mechanical frequencies are equal, see section V B 2. Even for instantaneous measurements, the agreement should still be better than one part in 2πk 0 N p (n + 1  2 ), which is a significant challenge.The most serious issue, however, comes from environmental decoherence.For both modulated and unmodulated couplings, we find that the relaxation rate must be bounded by τ R ≤ 2k B T Gm d 3 .This agrees with results previously reported in other scenarios, for example in two free oscillators [27], linearised optomechanics [28] and levitated nano systems [65], which supports the expectation that this is inherently a property of the noise model (57).That is, the constraint is independent of both the local dynamics and initial conditions -it cannot be improved upon with novel quantum control.To put this in context, even using the highest density materials (2.2 × 10 4 kg/m 3 ) and milliKelvin temperatures, the relaxation time needs to be longer than τ −1 R = Q/ω m ≈ 10 14 s.The requirement on the Q factor can be reduced with low frequency mirrors, however cooling then becomes more demanding.
There are two caveats here.In practice, equation ( 57) is likely to be too simplistic a description, and there is evidence that, at least in some optomechanical systems, the environment turns out to be non-Markovian [66].This is likely to have a further negative impact entanglement generation [67].On the other hand, there are also suggestions that the steady state limit of the full quantum Brownian motion model is more protected against loss of coherence [68].Whether this holds in the non-linear regime and for sufficiently short timescales -particularly in macroscopic systems where the discounted terms will typically be much smaller than those in (57) -remains to be seen.Alternatively, one may also look to weaker quantum correlations.However, in this case the motivation is less clear, as one no longer has the guarantee of monotonicity under LOCC (e.g.Discord [69]).
states with large variances for the coupling operators.For some states, the scaling in (A11) effectively extends to all values of the linear entropy.Setting Ôj = Nj , then one can readily confirm that for Fock state superpositions of the form |ψ j = 1 √ 2 (|0 j + |N j ), the linear entropy in (A1) is given exactly by, S L = 1 4 (1 − cos DN2 ) [71].Therefore ( 19) is satisfied with c = 1 4 cos −1 (1 − 4S L ), for 0 ≤ S L ≤ 1 2 , and in particular the maximum entanglement scales as D = π/N 2 .For fixed photon numbers, N p ≡ Nj = N/2, these states maximise the variance and can be considered optimal.This should come as no surprise -under the evolution (13) they lead to a superposition of two different displaced motions for the mechanical elements (i.e.oscillating cat states).The enhancements with N/2 from each cavity (equal to the standard deviations ∆ Nj ) are then roughly proportional to the maximum separation of the superposition (see also Example 1. in section VIII).
Note in this case, one can also readily calculate the entanglement through the von Neumann entropy of one of the reduced states, S = −Trρ 1 ln ρ 1 = − j λ j ln λ j , where λ j are the eigenvalues of ρ 1 .Under the assumptions above, it is straightforward to show, For a two level system, ρ 1 corresponds to a 2 × 2 matrix, and so the eigenvalues are given by the well known formula, where from (A12), As the entropy is maximised (but not necessarily maximal) when the separation between eigenvalues is largest, we again have the condition cos D(u In optomechanical systems, however, generating Fock state superpositions with large N is extremely difficult and a much higher photon number variance can be achieved by using coherent states.In this case, (19) does not extend all the way up to the point of maximum entanglement.This is already clear from the next term in (A10), which can be expanded to arbitrary order using the identity α| N n |α = B n (N p ), where B n (x) is a Bell polynomial.While this form can be valid for large N p (small D), its accuracy typically grows slowly for increasing S L .An alternative expansion can be found starting from the purity (A5) with |ψ 1 = e − Then setting X = D(p − q), with n = 0 in the identity (C9) (see below) leads to, which can alternatively be recast in the slightly more computationally convenient form, As long as N p is not too large, this allows a reasonably accurate evaluation of the linear entropy up to S L ∼ 1.We find that for most values, this turns out to be approximately a function of DN p only (assuming (∆ Nj ) 2 = N p is identical in each cavity).In practise, this means that we can still use (19) to estimate the entanglement time.A conservative choice for the measurable entanglement threshold is to assume S L = O (1), where now the proportionality constant needs to be determined numerically.For concreteness we choose the point corresponding to the maximum rate of increase the linear entropy, which is found to be at approximately |D(τ e )| = 1 . One then notes that in the macroscopic limit the second term in (B1) can be ignored, leading to This equation is Markovian, and could in principle have also been derived within the Lindblad framework.For two weakly interacting systems, each with their own bath, we can expect that the indirect coupling of one system to the other's environment will be very small.In the proposal considered here, this should hold for any realistic experiment, where γ g 1, and so it is reasonable to suppose a master equation of the form (57), where Υ j = 2m j k B T j γ R,j is the normal diffusion coefficient associated to each cavity, and for convenience we take Ĥren = Ĥ.This is justified as the mechanical frequency shifts in the expected regimes will typically be very small and therefore not greatly effect the entanglement dynamics.The extension from (B2) to (B3) has been validated both from first principles derivations, and numerically, for simple systems, such coupled harmonic oscillations [47], and shows good agreement in the small coupling limit.Furthermore, in the case of equal temperature baths, (B3) leads to essentially identical results to the strong coupling description, if positivity is enforced via the secular approximation.It should be noted that the same may not be true for the optical modes, where in particular we aim to work in the strong optomechanical coupling regime.In the present work, however, we will not address combined optical and mechanical noise.
For our purposes, it is sufficient to solve (B3) for the reduced cavity state alone.Here we make use of a straightforward generalisation of the approach in [48].For completeness, we repeat the relevant steps below.First, we transform to the the c.o.m. and stretch modes so that Ĥ is diagonalised.When the cavities are identical, i.e.Υ 1 = Υ 2 then (B3) can be rewritten as, where Υ = 2mk B T γ R .We then define the density matrix where |m, n = |m 1 |n 2 are field states of the two optical cavities.Thus the components of the reduced cavity state are given simply by, Using the projection (B5), we can re-write (B4) as, where we have used that Ĥ is diagonal in the photon number basis and we defined Ĥm,n = m, n| Ĥ|m, n as the effective Hamiltonians acting on the normal mode subspace, with, The next step is to move to the interaction picture by defining ρI D (t) = Û p,q † (t)ρ D (t) Û m,n (t), where the Û m,n (t) are defined through, and since Ĥ is photon number conserving, Û m,n (t) = m, n| Û (t)|m, n .Thus we have, We now multiply from the left and right by Û p,q (u) and Û m,n † (u), respectively and take the partial trace over the mirror modes, where we have used the cyclic property of the trace in the second line.The double commutators can be evaluated using where and The first set of phase factors do not contribute to xµ (t), and so we find, and so the double commutators become Substituting into (B12) we find, where, This can be solved by a simple integration, and so (setting u = t) we find, where ρ I D (0) = ρ D (0) = p, q| ρ(0) |m, n .This leads to an effective heating of the mechanical modes and a thermal contribution to reduced cavity state which no longer vanishes at the respective decoupling times of the optical and mechanical states.To see this explicitly we note from Appendix D that under unitary evolution, the components of the (reduced) cavity state when the mechanics is initially in a thermal state are given by, where For the modulation frequency Ω n = (1 − 1/n)(ω c + ω s ) and φ 1 = π/2, κ dec µ (t q ) can be evaluated to zeroth order in γ g as for µ = c, s.In section V A it was noted that the idealised limiting state of the final (entangled) optical fields will be of the form While many of the simplest continuous variable criteria fail to identify entanglement in states of this form, Shchukin and Vogel [40] introduced an infinite series of inequalities for which the violation of any is sufficient to witness entanglement.These can be conveniently written as the determinant of a matrix of moments of the state, where negativity implies the violation of one set of inequalities.The simplest which is suitable for our purposes is given by, For convenience we move to the Heisenberg picture, Ô(t) = Û † Ô Û , where in the symmetric cavity limit the evolution operator is given explicitly by ( 13) and ( 27), Together with D † µ (X Nj )â j D µ (X Nj ) = D µ (X)â j , and e −X N 2 âe X N 2 = e X(2 N +1) â we find, Where φ 0 = −ω 0 t + A + B and, are operators acting on the mechanical modes.Here we consider two scenarios: 1) The mechanical modes initially in a coherent state β µ , and 2) The mechanical modes initially in a thermal state ρ th,µ .For each case the expectation value of the displacement operator is given by, So we have where βµ = −(K Nµ β * µ − K * Nµ β µ ) and The final ingredient are the expectation values with respect to the cavity states, which we take to be coherent with parameters α 1 and α 2 for the two cavities, respectively.Using we have, These expressions are generally too complicated to gain much insight.However, by considering the form of the determinant (C1) we observe that the phase φ 0 ultimately vanishes together with the imaginary exponential terms in (C7) and (C8).Similarly, only the modulus of α 1 and α 2 is important and so from symmetry of the set-up this implies that the optimal choice of initial cavity states is α 1 = α 2 ≡ α and |α| 2 = N p .With these restrictions, the determinant witness (C1) can be written as, where κ th µ = |K Nµ | 2 (n µ + 1/2).We note, however, that the B terms should not in principle contribute to the entanglement, even though they have a non-trivial effect on the witness.An analogous situation arises in the estimation of the gravitational field by a single optomechanical cavity, where it is known that the quantum Fisher information (QFI) is independent of the corresponding B coefficient at the decoupling time, however certain POVMs only saturate the QFI when this is an integer multiple of 2π [21].This is occurs when using homodyne measurements, for example, and typically require k 0 ∈ Z. Interestingly, we find from a numerical analysis, that the optimal values of B(t q ) for the witness are integer multiples of π.Considering a measurement time t q = qnπ/(1/ω c + 1/ω s ), with the optomechanical coupling modulated at frequency Ω n = 1 2 (1 − 1/n)(ω c + ω s ), we find from B = B c + B s that, 2 .
(C12) Thus in the large n limit, we also see that approximately integer values of k 0 and are favourable [74] Adopting these B(t q ) values leads to immediately to the witness in the main text (28).The prefactor (generically N 3 p ) should not be attributed to the entanglement scaling estimated in section IV B. Instead this term acts like a signal gain, which should be offset against a proper noise analysis of measurements of the moments used in (C1).
When the initial mode temperatures are equal, we find that κ th c (t q ) = κ th s (t q ) to first order in γ g , κ th µ (t q ) ≈ (2πnγ g ) 2 k 2 0 (2n − 1) where in the second line we have assumed n large.Thus W 1 (t q ) = 4N 3 p e −4y e 2y sinh 2 (y) − z 2 , (C14) where y = 2N p z 2 + κ th (t q ) and z = sin(D(t q )/2).Taking the derivative with respect to z, the minimum will satisfy (34), (see (33)), which immediately leads to (38) under the substitutions N p → N p + Γ th 0 /2 and κ th → 0. However, for modulated coupling, where κ th = Γz 4/3 , finding a formal solution to the minimum condition is difficult, and so we instead resort to an expansion of (C15) under the assumption z min , N p z 2 min 1, and N p 1 (which are all consistent with the expected experimental regimes).In order to preserve the correct behaviour as κ th → 0, it is convenient to first multiply (C15) by (3/2)(4N p z min + κ | z min )/z min , before expanding the resulting expression to second order in z min .A slight rearrangement leads to, (C20) Transforming back to z min leads to (42) in the main text.This corresponds to an optimal entanglement time of, We now expand the exponential terms in W 1 (t 1 ) in powers of 1/N p , ignoring overall terms proportional to any inverse powers of N p (which is consistent with our final bound on Γ).This leads to, In the limit N p 1 the first term can be neglected, and so we again recover the asymptotic noiseless witness minimum of −N p /4.The solution for Γ such that (C22) is half the optimal value is then approximately, Γ 0.33N M mpnq = a m a * p b n b * q e −i(ω 0 t−A)(m+n−p−q) e iB(m 2 +n 2 −p 2 −q 2 ) e iD(mn−pq) e and nµ = 1/(e ω µ /k B T − 1) is the average thermal occupation of the transformed mechanical modes.Note that as entanglement is invariant under local operations, we can equivalently consider the state ρcav (t) defined by the coefficients, Mmpnq = a m a * p b n b * q e iD(mn−pq) e where the local actions on the optical states in ( 13) have been ignored.Further simplifications can also be made if we restrict our attention to short times around the optimal measurement time, t q − ∆t ≤ t ≤ t q + ∆t.This is consistent with our expectation of a small measurement window, and for nc = ns = n, an optimal strategy is to choose the modulation frequency Ω n = (1 − 1/n)(ω c + ω s )/2, with t q = qnπ(1/ω c + 1/ω s ).The difference between |K Ns (t)| to its 4 × 4 matrix representation (ρ T 2 cav (t)) (m,n),(p,q) = Mmpqn , where the indices run over (j, k) = (0, 0), (N, 0), (0, N ), (N, N ).The trace norm is then found from the sum of the eigenvalues of (ρ T 2 cav (t)) † (ρ T 2 cav (t)), leading to, This expression shows that the logarithmic negativity is non-vanishing only in regions where sinh(κN 2 ) < sin DN 2 2 .These define the available measurement windows in which entanglement can be detected by measurements on the optical fields.In particular, for small arguments of the sinh and sin functions, the condition takes the form κ ≤ D/2 which coincides with equation (32).A repeat of the analysis leading to (49) is straightforward, though in this case κ th max ≈ sinh −1 (1)/N 2 and so from (45), δt ≈ 2 sinh −1 (1) Note that as one can expect the optimal measurement time now scales as t q ∼ 1/N

2 FIG. 1 .
FIG.1.Two Fabry-Pérot cavities brought in close proximity.The moving elements interact via their gravitational fields.To second order this leads to an x1 x2 coupling and allows the generation entanglement between the mechanical modes.This is exchanged back and forth with the cavity fields, which themselves become entangled around multiples of the mechanical period 1/ω m .

Appendix C :
Evaluation of Determinant-based Entanglement Witnesses for Continuous Variable States

4z 2 min 1 2 2
+ e 2y min = 1 + 2z min 4N p z min + κ th | z min .(C15)For κ th | z min ≈ 0, a formal solution for z min can be written as, where W (x) is the Lambert-W function defined as the solution to W (x)e W (x) = x.For large x, we can approximate W (x) asymptotically asW (x) = ln x W (x) ≈ ln x ln x (C17)then to first order in 1/N p , the term W (e +N p +2κth N p ) = 1 2 + N p − 1 2N p [1 − 4(N p − 1)κ th ],and soz min ≈ 1 − 4(N p − 1)κ th the simple analysis provided in the main text.When κ | z min is significant, it is convenient to first write κ explicitly as a function of z.For = 0, κ th (z) ≈ Γ th 0 z