Optimal estimation of time-dependent gravitational fields with quantum optomechanical systems

We study the fundamental sensitivity that can be achieved with an ideal optomechanical system in the nonlinear regime for measurements of time-dependent gravitational fields. Using recently developed methods to solve the dynamics of a nonlinear optomechanical system with a time-dependent Hamiltonian, we compute the quantum Fisher information for linear displacements of the mechanical element due to gravity. We demonstrate that the sensitivity can not only be further enhanced by injecting squeezed states of the cavity field, but also by modulating the light--matter coupling of the optomechanical system. We specifically apply our results to the measurement of gravitational fields from small oscillating masses, where we show that, in principle, the gravitational field of an oscillating nano-gram mass can be detected based on experimental parameters that will likely be accessible in the near-term future. Finally, we identify the experimental parameter regime necessary for gravitational wave detection with a quantum optomechanical sensor.


I. INTRODUCTION
Precision measurements of gravitational effects allow for new technological advancements and for hitherto uncharted regimes of physics to be explored. In particular, the recent detection of gravitational waves by the Laser Interferometer Gravitational-Wave Observatory (LIGO) collaboration [1] has enabled the establishment of the field of gravitational astrophysics [2]. At the other end of the scale, fundamental tests of gravity using optomechanical systems have been proposed, including tests for gravitational decoherence [3,4], and measurements of the gravitational field from extremely small masses in quantum superpositions. Performing these experiments could help probe the overlap between quantum mechanics and the low-energy limit of quantum gravity [5][6][7][8][9][10][11]. Both endeavors are set to benefit from advances in quantum metrology [12], where the inclusion of non-classical states promises to push the sensitivity even further. This is already the case for LIGO, where the addition of squeezed light has significantly reduced the noise in the system [13].
Cavity optomechanics [14] represents a promising platform for developing high performance quantum sensors. These systems consist of light interacting with a small mechanical element, such as a moving-end mirror [15] or a levitated sphere [16]. In recent years, a diverse set of platforms including systems with Brillouin scat-tering [17][18][19], nanomechanical rotors [20], whisperinggallery-mode optomechanics [21,22] and superconducting devices [23] have been studied from both a theoretical and experimental point-of-view. When a mechanical mode is cooled to a sufficiently low temperature, it enters into a quantum state, which allows for properties such as entanglement and coherence to be used for the purpose of sensing [24].
Precision measurements of gravitational accelerationalso known as gravimetry -with quantum optomechanical systems in the nonlinear regime have been theoretically considered for measurements of constant gravitational accelerations [25,26]. However, constant signals are experimentally difficult to detect as they cannot be easily distinguished from a random noise floor. On the The influence of a time-dependent gravitational acceleration on a Fabry-Pérot moving-end mirror optomechanical system. A small source sphere with mass mS oscillates with frequency ωg and creates an oscillating gravitational field, which drives the center of mass motion of the mechanical part of the optomechanical system with frequency ωm. other hand, until recently it was not known how to solve the dynamics of fully time-dependent systems in the nonlinear regime. This prevented the careful analysis of the measurement of time-varying signals, which can provide significant advantages through the use of resonant effects.
The closed dynamics of time-dependent optomechanical systems was recently solved in [27], and a general expression for the sensitivity of an optomechanical system with a time-dependent coupling and time-dependent mechanical displacement terms was derived in [28]. In this work, we go beyond the results in [28] by deriving fundamental bounds to measurements of time-dependent gravitational fields and considering enhancements to the fundamental sensitivity. We apply our methods to three specific examples: generic gravimetry of oscillating fields, the detection of the gravitational field from small oscillating source masses, and the detection of gravitational waves. The computations are performed for both coherent and bright squeezed states of the light. We ask whether the intrinsic properties of the optomechanical probe system, such as the form of the light-matter coupling, can be employed to further enhance the sensitivity. This is motivated by the fact that the nonlinearities in the optomechanical coupling can be significantly enhanced by either separately or jointly modulating the mechanical frequency and the light-matter coupling [29,30]. Such modulations have been demonstrated e.g. in nanomechanical setups [31] or with levitated nanoparticles, such as in hybrid-Paul trap systems [32][33][34]. We find that such a modulation, performed at or close to resonance, significantly enhances the system sensitivity. A similar result holds when the trapping frequency is modulated at parametric resonance (twice the mechanical frequency).
To relate our scheme to realistic laboratory measurements, we also compute the sensitivity bounds for homodyne and heterodyne detection of the cavity state. While it was known that homodyne detection is optimal for constant gravitational fields and coherent states of the light in the cavity [25], here we show that it remains optimal for time-varying gravitational fields using initially coherent states of the optical mode (referred to as 'optics' for short in the following), as well as asymptotically optimal for squeezed states.
The work is structured as follows. In Section II, we introduce the optomechanical Hamiltonian and demonstrate how an external gravitational source enters into the dynamics. We outline the solution to the dynamics in Section II. Following that, we compute the quantum Fisher information for initial coherent states and squeezed states in Section III and discuss when the optical and mechanical degrees-of-freedom disentangle, since this allows us to focus exclusively on the sensitivity based on measurements of the cavity state. We then present our main results, which include expressions for the fundamental sensitivities in Section IV. Next, we consider homodyne and heterodyne measurement schemes in Section V. Finally, we apply our results and consider realistic parameters for three measurement schemes in Section VI: (i) generic gravimetry of time-dependent signals, (ii) detection of gravitational fields from small masses, (iii) and detection of gravitational waves. The paper is concluded with a discussion covering some of the practical implementations of an optomechanical sensor in Section VII and some closing remarks in Section VIII.

II. THE SYSTEM
The standard optomechanical Hamiltonian for a single interacting optical and mechanical mode is given bŷ whereâ,â † andb,b † are the creation and annihilation operators of the cavity field and mechanical oscillator, respectively, satisfying [â,â † ] = 1 and [b,b † ] = 1, where and ω c and ω m are the optical and mechanical frequencies. The light-matter coupling is denoted k, and its precise form depends on the experimental platform in question 1 . We consider the case where an external gravitational signal affects the mechanical element, which gives rise to an additional potential term in (1). However, this is not the only change to (1) that we consider, as will become clear later. By expanding the gravitational potential to first order, we obtain the familiar expression mg(t)x m , where g(t) is a time-dependent gravitational acceleration andx m = x 0 b † +b is a linear displacement of the mechanical element, with x 0 = /2mω m the zero-point fluctuation. While generic time-dependent signals can be explored using the methods in [28], here we restrict our analysis to gravitational signals g(t) that are sinusoidally modulated around a constant acceleration, which accounts for all three examples that we model in this work.
We write the gravitational acceleration as g(t) = −g 0 (a + cos(ω g t + φ g )) , where g 0 is the overall amplitude of the acceleration, φ g is an arbitrary phase, a is a dimensionless constant contribution, is a dimensionless oscillation amplitude, and ω g is the angular frequency of the signal. This allows us, for example, to model the gravitational field from an oscillating spherical source mass, as illustrated in Figure 1, where a = 1 and = 2δr 0 /r 0 , with δr 0 being the amplitude of the time-dependent oscillation and r 0 the mean separation (see the derivation in Appendix A). We can also use (2) to model gravitational waves (or a set-up mimicking their effects using, for example, moving masses [35]). To do so, we set a = 0 so that only the oscillating part of the gravitational acceleration remains.
It is well-known that resonances in physical systems can be used to further enhance certain dynamical effects. We therefore make a total of three changes to the standard optomechanical Hamiltonian (1): (i) We add a gravitational term g(t), (ii) we promote the standard constant optomechanical coupling k to a time-dependent one, and (iii) we let the mechanical frequency change as a function of time. The change of the mechanical frequency (iii) is achieved through the addition of a term D 2 (τ )(b † +b) 2 , which allows us to keep the mechanical mode operators fixed. The Hamiltonian in the frame rotating with the optical field then becomeŝ where we adopt a rescaled time parameter τ = ω m t, and the linear gravitational displacement term D 1 (τ ) becomes (given (2)) where Ω d1 = ω g /ω m , φ d1 = φ g , and where we now identify The optomechanical coupling k(τ ) depends on the specific system under consideration. For example, for a Fabry-Pérot cavity with a mechanical oscillator mirror, the coupling is a constant, k(τ ) ≡ k 0 given by k 0 = x 0 ω c /(Lω m ) [36], where L is the length of the cavity. For levitated dielectric spheres, the coupling takes the form k 0 = P k c x 0 ω c /(2ω m V c 0 ) [37], where P is the polarizability of the sphere, given by P = 3V 0 (ε − 1)/(ε + 2), with volume V , relative permittivity ε, and the cavity mode volume V c . Furthermore, 0 is the vacuum permittivity, and k c = 2π/λ c is the wave number of the light field. A modulated spring constant k(τ ) is experimentally feasible for Fabry-Pérot systems by positioning an electrode with a time-varying voltage close to the cantilever [38]. For a levitated nanosphere, a similar modulation arises from the natural micromotion that occurs for certain hybrid Paul-trap setups [32,39]. We later show that a modulation of the light-matter coupling can be used to enhance the sensitivity of the system for measurements of gravitational fields. Our goal now is to solve the dynamics generated by (3). The full solution was developed in [40] and [27]. We briefly review the results here. In general, the timeevolution operator is given by the time-ordered expo-nentialÛ (τ ) = ← − T exp −i τ 0 dτĤ(τ ) . By using an approach akin to transforming to the interaction picture, U (τ ) can be written as the product wherê Here,Û sq encodes both the free evolution of the mechanical subsystem as well as the term multiplied by D 2 (τ ), whileÛ NL contains the remaining nonlinear light-matter interaction term and the gravitational displacement term.
Next, we use a Lie algebra approach to write the remaining time-evolution operatorÛ NL (τ ) as a product of unitary operators. This method was first proposed by Wei and Norman in 1963 [41] and has since been used to solve the dynamics of a large variety of systems [42][43][44][45]. We identify the following Lie algebra of generators, which is closed under commutation: Similarly, it is possible to find a Lie algebra that generates the time-evolution encoded inÛ sq . It is made up of the following operators [27]:N b ,B Identifying the Lie algebra enables us to write down the following Ansätze for the two time-evolution operators [27] By now equating the two Ansätze (9) with their respective expressions in (7) and differentiating on both sides, we can use the linear independence of the operators to obtain a number of differential equations. Solving these, we find that the F coefficients are given by integrals shown in (B4) in Appendix B, and the J coefficients are similarly given by the expressions in (B8). For explicit expressions of the functions k(τ ), D 1 (τ ) and D 2 (τ ), it is then possible to solve the system either exactly or numerically.
In this work, we draw on analytic and perturbative solutions developed in Refs [27] and [28], which are briefly outlined in Appendix B.
A. Initial states of the system It is well-known that the fundamental sensitivity of a detector depends on the initial state of the system, and that significant enhancements can be gained through the use of non-classical states. For optomechanical systems, ground-state cooling has been demonstrated for a number of platforms [24,[46][47][48], however, the most realistic and practical state of the mechanical oscillator is a thermal state. The total initial state of the system iŝ where |ψ c is the initial optical state of the cavity and the parameter r T is defined by the relation r T = tanh −1 (exp[− ω m /(2 k B T )]), for which k B is Boltzmann's constant and T is the temperature.
In this work, we consider two different cavity states: (i) A coherent state |µ c (accessible through laser driving), whereâ |µ c = µ c |µ c . The average number of photons in the cavity is |µ c | 2 .
It is known that a Fock state superposition given by (|0 +|n )/ √ 2, whereâ †â |n = n |n can be used to maximise the sensitivity of the system for a given maximum excitation n [54,55]. However, it is difficult to prepare states with large n (currently, n = 4 has been experimentally demonstrated [56]), and we therefore focus on (squeezed) coherent states in this work.

III. QUANTUM METROLOGY OF LINEAR DISPLACEMENTS
We are interested in the fundamental limits that optomechanical systems can achieve when sensing displacements due to gravity. For this purpose, we turn to tools in quantum metrology.

A. Quantum Fisher information
In general, quantum metrology provides an ultimate bound on the precision of measurements of a classical parameter θ. If θ parametrises a unitary quantum chan-nelÛ θ , it is coded into the state asρ θ =Û θρinÛ † θ [12]. Then, given a specific input stateρ in , it is possible to compute the quantum Cramér-Rao bound (QCRB): Var(θ) ≥ (M I) −1 , where I is the quantum Fisher information (QFI) for the parameter θ and M is the number of measurements, or input probes [57]. The QCRB bound is optimised over all possible (POVM) measurements and data analysis schemes with unbiased estimators, and can be saturated in the limit of M → ∞.
In our case,Û θ is the time-evolution operator in (6), which results from a gravitational signal affecting the optomechanical system. The general form for the global QFI for the Hamiltonian (3) was computed in [28].
We are interested in estimating parameters that appear in the displacement function D 1 (τ ), which arises from the gravitational signal. We therefore pick d 1 in (4) as our fiducial estimation parameter, and by the chain-rule, we can choose to estimate any parameter that appears in d 1 . With this choice, only three dynamical coefficients, FN a , FB + , and FB − , contain the function D 1 (see (B4) in Appendix B), meaning that all other coefficients are zero when differentiated with respect to d 1 .
It follows from equation (9) in [28] that the operator H d1 is given bŷ where B and C ± are coefficients defined by The global QFI takes the form (see the derivation in Appendix C): where the variance ofN a , (∆N a ) 2 ≡ Var(N a ) = N 2 a − N a 2 , and where the bracket · denotes the expectation value with respect to one of the initial states presented in Section II A. For the coherent state and the squeezed state, we find (see Appendix C): where we recall that r and ϕ are the squeezing parameters given by ζ = re iϕ . The angle ϕ is defined with respect to the coherent state phase. The case of coherent states (r = 0) was considered previously in [25,26,28]. For coherent states, a higher photon number |µ c | 2 yields a better sensitivity. For squeezed coherent states, the QFI is maximised when e iϕ 2 µ c is purely imaginary, and when the photon number |µ c | 2 and r are maximised. In each case, the increase in sensitivity is not without cost, as there are certain restrictions to how much the mechanical element can be displaced. See Section VII B for a discussion of these restrictions, where we also propose order-of-magnitude limitations for the parameters of the cavity field.
Once the QFI has been computed, we can obtain the optimal measurement sensitivity through the QCRB. Given the dimensionless expression (4), we use the chainrule to find that the sensitivity ∆g 0 to the gravitational amplitude g 0 (see the expression in (2)) is In Section V, we consider the classical Fisher information (CFI), which provides the sensitivity given a specific measurement. We now turn to the question of optimal timing of the measurement.

B. Disentangling of the optics and mechanics
For Hamiltonians such as (3), it is well-known that the optical and mechanical subsystems evolve into an entangled state [60], however, for particular choices of the dynamics, we find that there are certain times when the two systems end up in a separable state. This is a consequence of the unitary dynamics, and we refer to these times as τ sep .
In an experiment, it is often the case that an observable on the cavity state is measured. If we can identify the disentangling conditions and hence τ sep , we can immediately compute the QFI of the cavity state at these separation times. It was also previously found that the global QFI peaks when the states are separable, and that the noise contained in an initially thermal mechanical state also does not affect the sensitivity at this time [25,28].
From (6), we note that only the exponentials with FN aB+ and FN aB− mediate entanglement between the cavity field and the mechanics, since their accompanying generatorsN aB+ andN aB− encode an interaction between the light and mechanical oscillator (referred to as 'mechanics' for short in the following). We therefore construct the function KN a = FN aB− +iFN aB+ , and express a sufficient condition for separability as When this condition is fulfilled, the full time-evolution operatorÛ (τ ) factorises into terms that act exclusively on the optical and mechanical states. The part acting on the cavity state is given byÛ c (τ ) = For later, we note that, when applied to a coherent state |µ c , the last exponential induces a phase, such that the new coherent state parameter isμ c = µ c e −iFN a . This definition will become useful to us when we discuss homodyne measurements of the cavity field in Section V.
The advantage of identifying the conditions for |KN a | 2 = 0 is that we can derive an analytic expression for the fundamental sensitivity that can be achieved by measuring the cavity state. We also do not have to concern ourselves with any contributions from the initial thermal mechanical state. The QFI of the optical state is then simply (from (13) and (14)): where we use the subscript 'c' to denote the fact that this refers to the QFI of the cavity state only.
To determine when the condition in (18) is satisfied, we must evaluate the expression for a given dynamics. Firstly, we note that the form of the gravitational acceleration (determined by the function D 1 (τ )) does not affect the entanglement between the systems. This is because D 1 (τ ) does not feature in the integrals for FN aB+ and FN aB− (see the expressions in (B4)). In contrast, the optomechanical coupling k(τ ) and the squeezing function D 2 (τ ) completely determine the times τ sep at which the two systems become separable. For a constant optomechanical coupling k(τ ) ≡ k 0 , the states reach their maximum entanglement at τ = π, after which they return to a separable state at τ sep = 2π [25,60]. We can prove this explicitly by computing FN aB+ and FN aB− for a constant coupling, and we find that |KN a | 2 = 2k 2 0 (1 − cos(τ )), which vanishes when τ is a multiple of 2π.
When the coupling k(τ ) is time-dependent, however, the behaviour of the system -and the entanglementbecomes richer. As we are interested in whether a modulated coupling can lead to resonance type enhancements, a natural choice is to assume it takes on the form [30] where Ω k is the oscillation frequency divided by ω m and φ k is an arbitrary phase. For zero mechanical squeezing (D 2 = 0), the F -coefficients are given in (B10), and KN a is given in (D1). When the optomechanical coupling is modulated at resonance with Ω k = 1, we find that the light and mechanical oscillator never disentangle. This means that we cannot ignore the mechanical contribution to the QFI, and since computing the QFI for a reduced state is challenging, we resort to the global expression in (14) as an upper bound.
A key observation however, is that for specific choices of the coupling modulation frequencies, the light and mechanics do disentangle at certain points in the evolution. In Appendix D we prove that when the frequencies take ) for a modulated optomechanical coupling (specified in (20)) with different frequencies Ω frac = 1 + 2n1/s and with k0 = n1 = 1. The resonant case Ω k = 1 (not shown) never evolves into a separable state, while |KN a | 2 vanishes at multiples of τ = sπ for Ω frac = 3 (blue solid line), Ω frac = 2 (green dashed line), and Ω frac = 5/3 (magenta dotted line). (b) is a plot of 2 log |KN a | for a modulated squeezing with d2 = 0.01 and a constant optomechanical coupling k0 = 1. At no point within the shown time interval does the states evolve into a separable state. on a fractional form, Ω frac = 1 + 2n 1 /s , for n 1 and s integers (s positive), the subsystems decouple at times that are multiples of τ sep = sπ. This means that the QFI for the cavity state is given again by (19).
Finally, for a mechanical frequency modulated with D 2 = d 2 cos(2τ + φ d2 ), we find no point where the system is completely separable (see Figure 2b).

IV. GRAVIMETRY OF TIME-DEPENDENT GRAVITATIONAL FIELDS
We are now in a position to evaluate the QFI explicitly for a number of cases of interest. Throughout, we assume that the gravitational signal is given by the timedependent expression in (4). Further, we keep the optomechanical coupling constant for now with k(τ ) ≡ k 0 , and we assume that the mechanical squeezing is zero. In Section III B, we showed that for this choice of dynamics, the light and mechanical oscillator disentangle whenever the time τ is a multiple of 2π.
We therefore find that (see Appendix C 2), at resonance with Ω d1 = 1, and at time τ = 2nπ with n integer, the global QFI becomes which is maximized over φ d1 for φ d1 = π. This is a phase relation between the driving signal, which excites oscillations of the mechanics, and the light-matter coupling term, which fixes the decoupling times. The much longer form of the QFI for a general frequency Ω d1 is shown in (C20) in Appendix C.
In a classical setting, we expect that driving the mechanical element on resonance will rapidly increase its oscillation amplitude, which means that it becomes easier and easier to detect its displacement. We do see this increase in the ∼ n 2 -scaling of the second term of (21). However, this term is usually small compared with the first term, since both scale with ∼ n 2 and the first term scales with the photon number variance (∆N a ) 2 .
To focus on this point, we consider the cavity state QFI (19) at times when the light and mechanics evolve into a separable state. For a purely oscillating field with a = 0, the local QFI for measurements of the cavity field becomes When = 1, this is in fact smaller than the constant driving scenario (with a = 1) by a factor of 4. So, while resonant driving does increase the global and local QFI over time, as one would intuitively expect, this is primarily through the amplitude change of the mechanical element. As such, it does not translate directly to observations on the cavity field. It turns out, however, that an analogous enhancement can be passed to the field provided modulations are introduced to the system in a different way. In this section we consider two additional methods by which this can be done: Modulating the optomechanical light-matter interaction, and modulating the trapping frequency. We return to (22) later on and use it to compare the effects of the enhancements.

A. Enhanced sensing through optomechanical modulation
We are interested in whether the form of the lightmatter coupling k(τ ) can be used to enhance sensitivity of the system. Such a time-modulated coupling has been experimentally demonstrated [38,61,62]. We specifically consider a coupling of the form shown in (20). The global QFI for a resonant gravitational signal at arbitrary times can be found in (C21), and it is dominated by terms proportional to n 4 for large τ when (φ d1 −φ k )/π is not an integer. Should coherence be maintained for long periods of time, the resonantly modulated coupling leads to rapid increases in the measurement precision.
For mechanical resonance (Ω k = Ω d1 = 1), we noted before that the light and the mechanics do not disentangle at all. This means that the QFI is global at all times, and therefore does not necessarily reflect the sensitivity that could be realistically obtained in the laboratory through measurements of the optical state. For multiples of the mechanical period, τ = 2nπ, the global QFI becomes, The full expression for arbitrary τ is given in (C22). We note that the term multiplied by provides an additional scaling with n 2 , leading to an overall scaling of n 4 . Such an enhancement is only present when the gravitational field is oscillating with nonzero , and indicates that the two resonances (the gravitational signal and the optomechanical coupling, resonant with the mechanical frequency) constructively enhance the sensitivity. This is optimised when φ d1 − φ k = π/2. Furthermore, the term multiplied by a is maximised for φ k = 0, so we can choose φ d1 = π/2 to optimise the expression. For a purely oscillating gravitational field a = 0 and a large temperature r T → ∞, then setting φ d1 = π/2 and φ k = 0, simplifies the expression in (23) to which, compared with a constant coupling, is an improvement of ∼ n 2 π 2 /4 for purely oscillating fields (22). The global QFI is generally not accessible in an experimental setting, since it is difficult to measure the mechanical element directly. However, we saw in Section III B that the light and mechanical oscillator become separable for very specific choices of the frequency Ω k , which we referred to as the fractional frequencies Ω frac . With this choice, we compute the local QFI for the cavity state with Ω k = Ω d1 = Ω frac 2 . Using the expression for the cavity state QFI in (19), we find that when Ω frac = 1 + 2n 1 /s, with s > being a positive integer and integer n 1 = 0 with 2n 1 /s > −1, the QFI becomes, at τ sep = q sπ, where q is a positive integer: For a purely oscillating signal with a = 0, the optimal choice of phases φ d1 − φ k = 0, and n 1 = −1 (which means that Ω frac = 1 − 2/s and s ≥ 3), we find that the optimal choice for q and s for a given disentangling time τ sep = q sπ is to maximise s which implies q = 1. Then, the QFI becomes Equation (26) is one of the main results in this paper, since it provides a sensitivity that can be realistically achieved from measurements on the cavity state alone. To see how well the enhancement compares, we contrast I (Ω d1,k =Ω frac ) with (22). Note that n is the parameter giving the number of mechanical periods. The meaning of s is different; it is the parameter defining the fractional frequency. Using (26), and assuming that s is even, such that s = 2n, we find an improvement of ∼ n 2 /4 for s 1. For arbitrary times, we refer to Figure 3, which shows the general behaviour of the global QFI. The plot in 3a compares resonant gravimetry I (Ω d1 =1) with the enhancements I (Ω d1,k =1) and I (Ω d1,k =Ω frac ) obtained by including a time-dependent coupling k(τ ) for purely oscillating gravitational fields, and the plot in 3b shows the same quantities for a gravitational field with constant and oscillating parts. In both plots, we consider large temperatures with r T → ∞, (which minimises any additional information which can be gained from the mechanics), and set k 0 = (∆N a ) 2 = 1, since these are merely multiplicative factors in the QFI. We also choose the optimal phases for each setting, which are φ d1 = π for I (Ω d1 =1) , φ d1 = 0 and φ k = π/2 for I (Ω d1,k =1) , and finally φ d1 = φ k = π/2 for I (Ω frac ) .

B. Enhanced sensing through modulated mechanical frequency
The second enhancement we consider (separately from the above) is the inclusion of a mechanical squeezing term D 2 (τ ) b †2 +b 2 . We assume that it is periodically modulated with where d 2 is the amplitude, Ω d2 is the rescaled modulation frequency and φ d2 a phase factor. FIG. 3. The quantum Fisher information for detecting linear displacements with a modulated optomechanical system. We set rT → ∞, k0 = 1, µc = 1, and re iϕ = 1 in both plots. (a) shows the QFI for a purely oscillating gravitational field with a = 0 and = 1. We compare the global QFI for a resonant gravitational signal I (Ω d1 =1) in (21) (green solid line) with the enhanced global QFI for a modulation of the optomechanical coupling at resonance (Ω d1 = Ω k = 1) denoted by I (Ω d1,k =1) in (C22) (blue dashed line), and the enhanced QFI for fractional frequencies Ω d1 = Ω k = Ω frac = 1 + 2n1/s denoted by I (Ω frac ) in (C23) (dotted purple line), where we set n1 = −1 and s = 8 for this plot. The phases have been chosen such that they optimise the QFI for each case (see the main text). The QFI for a resonant coupling shows the strongest increase for later times, but the states never disentangle, which means that we can only upper bound the sensitivity for a measurement of the optical field alone. (b) shows the global QFI for a constant plus oscillating gravitational field with a = 1 and = 0.1, where we estimate the overall amplitude d1. The fractional frequencies no longer perform well because a, however the scaling with the parameter s can only be appreciated when comparing the curves for different fractional frequencies. Resonant gravimetry, denoted by I (Ω d1 =1) , increases smoothly but it is outperformed by modulated resonant gravimetry I (Ω d1,k =1) for large τ .
A term of this form can be generated by, for example, modulating the spring constant [30] or the trapping frequency of a levitated system [33,34]. In particular, in the levitated systems presented in [32][33][34], modulations of the light-matter coupling are always accompanied by a modulation of the mechanical frequency.
When Ω d2 = 2, this corresponds to a parametric amplification of the mechanical oscillation and leads to a squeezed state of the mechanics (see [63] for how this can be implemented experimentally). The perturbative solutions of the dynamics were found in [27], and are valid for d 2 1 and d 2 τ of order (at most) one. This means that we can only consider small values of d 2 , especially if we are interested in large times τ .
When the mechanical trapping frequency is modulated sinusoidally, the light and mechanics never disentangle, and we are therefore unable to consider the QFI of the optical state separately (see Section III B and Figure 2b). We therefore resort to the global QFI in (14). The modulation of the mechanical frequency leads to an enhancement of the QFI depending on the phases φ d1 and φ d2 , however the full expression is long and cumbersome. We refer to Appendix C 4, and instead find the optimal phase choice numerically. From the QFI plotted in Figure 4, we see that the choice of φ d2 = −π/2 and φ d1 = 0 maximises the QFI.
With this choice of phases, taking into account that d 2 τ ∼ 1 and d 2 1, the dominating term in the QFI is Compared with the QFI for resonant gravimetry without any enhancements in (22), the modulated mechanical frequency brings an improvement of ∼ (e − 1) 2 ∼ 3 when d 2 τ ∼ 1, and τ = 2π. This means that the addition of a modulated squeezing term can increase the sensitivity, but we are limited by our perturbative method in predicting its efficiency. The inclusion of a constant squeezing term D 2 ≡ d 2 is equivalent to changing the mechanical frequency as ω m → ω m √ 1 + 4D 2 . Since the dimensionful QFI scales with ω −5 m 3 , larger d 2 means that the QFI decreases.

V. HOMODYNE AND HETERODYNE METROLOGY OF LINEAR DISPLACEMENTS
While the QFI and the QCRB provide the ultimate limits to how well a parameter can be estimated, it is not immediately clear which measurements actually saturate this bound. Experimentally, one would almost always measure the optical state using a homodyne measurement, a heterodyne measurement, or photon counting.
The cavity field as present in our description is not directly experimentally accessible, although the contrary 3 The dimensionful QFI is proportional to k 2 0 , which in turn is proportional to ω −3 m . Furthermore, another factor of ω 2 m appears from the dimensionful factor given from the sensitivity in (17), which appears as a multiplicative factor in front of the QFI. The plot shows the exact QFI I (Ω d2 =2) for the sensing of a purely oscillating gravitational field (i.e. a = 0) for different values for the initial phase parameters of frequency modulation and gravity oscillation (φ d2 , φ d1 ). The gravitational field is modulated on resonance (Ω d1 = 1) and the mechanical frequency is modulated on parametric resonance (Ω d2 = 2). For the plots, we used the parameter values k0 = 1, = 1, d2 = 0.02, µc = 1 and no squeezing of the cavity field. However, k0, and µc appear only if rT is very large, as we assumed. Therefore, the only relevant parameter is d2, which defines the time-scale on which the effect of parametric driving becomes pronounced.
is commonly assumed in the literature. To build on these results, one would have to consider output fields leaking from the cavity, which we leave to future work. Instead, here we compute the classical Fisher information (CFI) for these ideal measurements on the cavity state, focusing on when the light and mechanics have disentangled (see Section III B).
When the light mode and mechanical oscillator are in a separable state, the local QFI generator reduces tô H d1 = BN a , where B is defined in (13). The optimal bound is given by the QFI in (19), and our aim is to investigate whether a homodyne or heterodyne measurement satisfies this.
The general expression for the CFI is where we henceforth denote all CFI quantities by I, rather than I, which we reserve for the QFI, and where p(x, d 1 ) = Tr ρ d1Πx is a probability distribution resulting from a measurement with a POVM elementΠ x . Assuming that the initial cavity-field state is pure (which in the settings we consider here is always true when the optics and mechanics are separable), we define the state a acts on the cavity state. Then, noting that the probability is given by p(x) = | ψ τ |x | 2 and 1 = dx |x x|, the CFI can be written, where The first term in (30) is relatively straightforward to calculate, however, it is generally difficult to perform the integral in (31). A particular simplification exists when x|N a |ψ τ is proportional to x|ψ τ . This occurs, for example, when the state at τ = τ sep is a coherent state, which can be guaranteed by choosing parameters such that the coefficient FN 2 a is a multiple of 2π at the disentangling time (see Appendix E for details). For mathematical convenience we will make this assumption in the remainder of this section, however it turns out that this special case is still sufficient to saturate the QFI for practical measurement schemes, unless the initial cavity state is squeezed (in which case the CFI still approaches the QFI for large photon number).

A. Homodyne measurements
We start by investigating homodyne measurements for coherent and squeezed coherent optical states, since these are standard measurements that are routinely performed in the laboratory.
In [25], it was shown that the QFI is saturated at τ = 2π by a homodyne measurement when the rescaled optomechanical coupling takes an integer value and when the gravitational acceleration is constant. The question is whether the homodyne measurement is still optimal when the gravitational field is time-dependent, and when a modulation of the optomechanical coupling is included.
In general, a homodyne measurement involves a measurement of the optical quadrature. The relevant POVM is given by |x λ x λ | where the state, |x = |x λ , is defined as the eigenstate of the operator, For an initial coherent state in the cavity, we show in Appendix E that the CFI is given by where B was defined in (13) (and thus contains the effects of modulating the coupling), and whereμ c = e −iFN a µ c . For matching choices of λ and µ c , the optimal value can always be found. When R[μ c e −iλ ] = 0, we find which coincides with the local QFI (19) for the cavity state. Therefore, we conclude that the CFI for homodyne measurements saturates the QCRB, provided that the phase λ can be optimally controlled.
A similar analysis can be performed when the initial optical state is squeezed. Adopting the convention |µ c , ζ =Ŝ(ζ) |µ c , where the squeezing parameter is given by ζ = re iϕ , we show in Appendix E that the maximum CFI (for a large photon number |µ c | 2 , such that it dominates over the vacuum contribution, and given the specific conditions in (E20)), is given by This is less than the maximum QFI (see the expression in (16)) by only a vacuum contribution. However, the CFI asymptotically approaches the QFI for large |µ c | 2 .
In general, however, the Fisher information can still be non-zero when µ c = 0. Here, we find the vacuum contribution a , and the F -coefficients are all evaluated at the time of separability. Similar to the QFI, the CFI reaches a maximum of I (hom) ζ = 2B 2 sinh 2 (2r) forφ − 2λ = ±2[nπ ± tan −1 (e −2r )]. However, for all but very small photon number (and large r) the optimal CFI is given by (35).

B. Heterodyne measurements
The heterodyne measurement case is somewhat more straightforward since the probabilities are calculated with respect to coherent states [64]. Replacing |x = |β , where |β is a coherent state, we find for |ψ 0 = |µ c the overlap appearing in R to be and so, The CFI for a heterodyne measurement is then 4 , 4 This corrects an erroneous factor of √ π/2 in [26].
which is half of the QFI (19) associated with the light field. For initially squeezed states, we have (see Appendix E 2), Similarly to the QFI, we find that when e − iϕ 2 µ c is purely imaginary, the CFI is maximised. However, it does not coincide with the QFI.

VI. IDEAL SENSITIVITIES FOR OPTOMECHANICAL SYSTEMS
In this Section, we use our results to obtain an orderof-magnitude estimate for the ideal sensitivity of gravity measurements. The sensitivities we derive below are merely indicative of the final sensitivities that can be achieved. We then briefly discuss squeezing of the cavity field and proceed to compute the fundamental sensitivity for three applications: generic accelerometry, sensing gravitational signals from small source masses, and detecting gravitational waves.

A. Fundamental sensitivities
We identify two key formulas from our results that provide the strongest sensitivities for the detection of time-dependent gravitational fields. Crucially, we limit ourselves to presenting sensitivities that we know can be achieved by homodyne measurements in the laboratory. This requirement rules out the enhancement that can be achieved when the optomechanical coupling is modulated at resonance and modulations of the mechanical frequency, simply because the system never evolves into a separable state. With our current tools, it is difficult to predict the sensitivity of a classical measurement on a mixed state, however this does not mean that high sensitivities cannot be achieved. We leave it to future work to explicitly explore those settings.
For measurements of the cavity state at multiples of τ = 2πn, the QFI for gravimetry of resonant gravitational fields in (21) leads to the sensitivity where we recall that m is the optomechanical mass, ω m is the mechanical oscillation frequency, k 0 is the optomechanical coupling, a is a constant contribution from the field and is the oscillation amplitude. We then allow the mechanical frequency of the optomechanical system and the optomechanical coupling to be modulated at the fractional frequencies Ω frac , which we identified in Section III B. We use the QFI expression in (26) to predict the following sensitivity for a measurement at τ = πs (s ≥ 3 being a positive integer), at which point the light and mechanical element are found to be in a separable state: where we have set = 1 and where we explicitly set a = 0. For bright squeezed states of the cavity field, ∆N a is maximised when µ c e iϕ/2 is fully imaginary, which can be achieved by assuming that µ c ∈ R and that ϕ = π/2. With this condition, we find that N a µc,ζ = µ 2 c cosh(2r) + sinh 2 (r), and (∆N a ) 2 µc,ζ = e 4r µ 2 c + sinh 2 (2r)/2. As mentioned earlier, it is common to report the squeezing in terms of decibel in experiments, which we call S dB . The relation between this quantity and r reads r = S dB /(20 log 10 e) [65]. Schemes for obtaining S dB = 10 have been proposed [66], which corresponds to r = 1.73. While the CFI for homodyne detection with squeezed states does not saturate the QCRB, it does so asymptotically as |µ c | 1 and small r.

Measuring oscillating gravitational fields
As the simplest application, we consider measurements of the oscillating part of a gravitational field. The constant part of the field (if present), can be absorbed into the system dynamics by letting the constant displacement of the mechanics be part of the initial state. This is equivalent to saying that we are performing a relative measurement of the gravitational field, where only the time-dependent part contributes. According to the equivalence principle, the sensitivity we derive here also applies to accelerometry measurements, when the optomechanical system is shaken with fixed frequency. Using the parameter values listed in Table I and considering a modulated optomechanical coupling, we find that the single-shot sensitivity predicted by equation (42) for measuring oscillating gravitational acceleration is ∆g 0 ∼ 1.4 × 10 −11 ms −2 . 5

Measuring gravitational fields from small oscillating masses
The interest in detecting gravitational fields from increasingly small masses stems from the desire to explore   (42) with Ω frac = 1 − 2/s = 9/10, where s = 20. For the chosen parameters, resonant gravimetry without modulation represented by equation (41) leads to a bound that is larger by about a factor 5 in comparison to the bound given by equation (42).
the low-energy limit of quantum gravity. If the gravitational field from superposed masses can be detected, it may, for example, be possible to examine how gravity behaves on these small scales [3,4,11,67]. An explicit setup for measurements of a miligram mass was proposed in [68]. We compute the fundamental bound for sensing gravitational fields from small source masses, which then allows us to place a limit on the masses that these systems can detect (for realistic source-detector separations). We refer to the expression for the gravitational potential in (A2) in Appendix A, where we have expanded the gravitational potential that results from small, timedependent perturbations from a moving spherical source mass. The resulting gravitational field oscillates around a constant value where g 0 g 0 . If the constant contribution can be measured, the most practical strategy would be to forgo any modulations of the coupling and consider the sensitivity given by (41). However, more realistically, it may lead to higher precision to estimate only the oscillating part (see for example [68]). In this case, the light-matter coupling can be modulated for an enhancement, and we use the expression in (42).
Given the values in Table I and a number of measurements M = 10 4 , we find that the maximum sensitivity for measuring the oscillating part of the gravitational field of a moving mass that can be achieved is ∆g 0 = 1.4 × 10 −13 ms −2 . For a spherical source mass oscillating with amplitude δr 0 at an average distance r 0 from the source such that the time-dependent distance is r(τ ) = r 0 − δr 0 cos(Ω d1 τ ), we find that (see Appendix A) the oscillating contribution to the acceleration is ≈ 2δr 0 Gm S /r 3 0 , where G is Newton's gravitational constant. We can solve for m S , and assuming that 2δr 0 /r 0 = 0.1, we find m S ∼ 200 ng given a dis-tance of 100 µm between the probe and source mass. At this distance, we expect the Casimir effect between the probe and source sphere to become noticeable, but this can potentially be remedied by shielding the system. We discuss this in Section VII D below.

Gravitational wave detection
Recent years have seen a surge in interest regarding the measurement of gravitational waves with novel setups, including proposals for detectors with superfluid Helium [69], Bose-Einstein condensates [70], and even interferometry with mesoscopic objects [71]. Here we investigate the feasibility of gravitational wave detection with an ideal cavity optomechanical system. Our approach is essentially the quantum analogue of the classical scheme presented in [72].
Compared with the previous section, here we focus on identifying the experimental parameter regimes needed to detect gravitational waves. The gravity gradient induced by a gravitational wave is given as G = ω 2 m h/2, and this is by far the dominant effect induced by a gravitational wave for optical resonator systems (see [35,73,74] for details of how deformable optical resonators can be described in a relativistic framework and how relativistic and Newtonian effects can be compared). Then, the differential acceleration between the two ends of the cavity system becomes g 0 = Lω 2 m h/2, and the error bound for gravitational wave strain h is given as Considering a single detector of 10 m length with the parameters given in Table II, we obtain ∆h ∼ 1 × 10 −21 .
Strains of the order of 10 −21 are expected for compact binary inspirals in the frequency range we considered here (see figure A1 of [75]). The time scale for a single measurement is τ /ω m ∼ 6 s, which is sufficiently short for several integrated measurements before the source leaves the considered frequency range ∼ 2 Hz (see figure 1 of [76]). The same sensitivity can be achieved by 10 4 sensors of length L = 10 cm (provided that Table II can be maintained).

VII. DISCUSSION
There are many practical aspects to building an optomechanical gravimeter, many of which are beyond the scope of this work. Features such as optical and mechanical noise are ever-present in experiments, and we briefly discuss these and other systematics in this section.   (42) with Ω frac = 1 + 2n1/s = 9/10. Similar numbers can be obtained considering a three-dimensional array of 10 4 detectors with cavity length 10 cm and 10 independent measurements with each detector to obtain for the total number of independent measurements M = 10 5 .

A. System parameters
In Tables I and II, we used example parameters to compute the ideal sensitivities from our results. We here discuss the feasibility of these parameters, and we identify a few features of different experimental platforms that appear beneficial for displacement sensing. Our aim is to provide a brief discussion of this topic rather than a comprehensive overview of the advantages and disadvantages for each experimental platform.
From our results, we see that a low mechanical frequency is beneficial for sensing. Considering the fact that k 0 ∝ ω −3/2 m , we see from (41) and (42) that ∆g 0 ∝ ω 3 m . At the time of writing, there are not yet many optomechanical experiments that have achieved ground-state cooling, and thereby operate in the quantum regime. Those that do (see e.g. [24,[46][47][48]) require high mechanical frequencies, which is therefore detrimental to sensing as envisioned here.
We do however identify a few platforms that lend themselves well for sensing of the kind explored in this work, although additional experimental progress is needed before these systems can operate in the nonlinear quantum regime. Crucially, levitated systems can achieve extraordinarily low mechanical frequencies; for example, particles levitated in a magnetic trap [77,78] can potentially reach mechanical frequencies as low as ∼ 2π×50 rad s −1 [79]. This type of system has in fact already been considered in the context of measuring constant gravitational acceleration [80]. A cavity could potentially be added to the magnetically levitated systems described in [78], which would allow the mechanical element to couple to the cavity field via a standard light-matter coupling of the form considered here. We note however that a lower frequency also requires the system to stay coherent for longer, which is of course challenging. Similarly low frequencies have been achieved with optically trapped nanoparticles [81].
Fabry-Pérot moving-end mirrors and membrane-inthe-middle configurations generally operate at higher mechanical frequencies (see e.g. [82]). However, we note that many of the features of optomechanical systems are interlinked, such as the mechanical frequency and the coupling constant. The sensitivities derived will therefore, in principle, be different for each unique setup.

B.
Restriction of the cavity-field parameters Our results remain valid as long as the dynamics of the system is well-approximated by the Hamiltonian (3). The standard optomechanical Hamiltonian is derived by assuming that the perturbation of the oscillator is small compared with a specific length scale of the system [36,83,84]. For a Fabry-Pérot cavity, the perturbation must be much smaller than the cavity length L [84,85], and for a levitated system, it must be smaller than or equal to the wavelength of the cavity light mode. This ensures that the radiation pressure force remains approximately constant [86].
If the mechanical oscillator is strongly displaced such that additional anharmonicities appear 6 in the Hamiltonian, our results can no longer be used to accurately predict the ideal sensitivity (that is however not to say that the system would perform badly as a sensor). We note that the system we consider is closed, and that the initial state corresponds to immediate radiation pressure on the mechanical oscillator. This effect is by far the largest contributing factor to the displacement (especially in the context considered in this work, where the gravitational effects are generally weak). To ensure that the oscillator is not displaced beyond the point at which the dynamics changes, we must consider restrictions to the parameters that determine the cavity field.
We introduce a generic length-scale l beyond which the extended Hamiltonian (3) is no longer valid. The nature of l will differ for each setup. Because the system is quantum-mechanical, we consider the probability of detecting the centre-of-mass of the mechanical element a certain distance away from the origin. This is wellcaptured by the expectation value x m and the standard deviation ∆x m , and we therefore require that they remain much smaller than the length-scale l at all times.
We derive explicit expressions for x m and ∆x m in Appendix F. The position of the mechanical oscillator is given byx m = x 0 b † +b , where we chose the equilibrium 6 Certain dynamics can still be solved, see e.g. [87]. position as the origin. When prepared in the groundstate, the general expression for the mean displacement and its variance as a function of time τ are given by where α(τ ) and β(τ ) are Bogoliubov coefficients that arise from the mechanical subsystem evolution, given in (B6), and where ∆(τ ) and Γ(τ ) are given by (F2) in Appendix F. With these expressions, we can identify the appropriate restrictions on N a and ∆N a for resonant gravimetry and the enhancement schemes presented Section IV. We here comment on the restriction for each scheme considered in Section IV: (i) For gravimetry without enhancements, we find that x m (τ ) oscillates with an amplitude 2x 0 (k 0 N a + d 1 a) about a mean displacement of the same size.
The restriction on the photon number is given as N a l/(2x 0 k 0 ). Analogously, the photon number standard deviation is restricted to ∆N a l/(2x 0 k 0 ).
(ii) If the light-matter coupling is modulated as a function of time, k(τ ) = k 0 cos(Ω k τ ), at resonance Ω k = Ω d1 = 1, the dominant term in x m is given by x 0 k 0 N a τ . Thus the condition on the photon number becomes N a l/(x 0 k 0 τ ). Analogously, the standard deviation is limited by ∆N a l/(x 0 k 0 τ ). This means that we are additionally limited by the integration time. The main difference to (i) is that both x m and ∆x m increase with time, which implies that the bound strengthens with τ .
(iii) Next, for Ω d1 = Ω k = Ω frac , and taking into account that τ sep = sπ, we find the condition N a , ∆N a πl/(x 0 k 0 τ sep ), which is larger by a factor π compared with the resonant case.
(vi) For modulations of the mechanical frequency with a constant coupling k = k 0 , we find the restriction N a , ∆N a l/(2x 0 k 0 (1 + e d2τ )), which holds for d 2 τ ≤ 1.
We note that the effect of the cavity field on x m can be canceled either by preparing the mechanical state in an appropriate coherent state, or by introducing an additional external potential that cancels the effect of the radiation pressure. When the light-matter coupling k(τ ) is modulated (which enhances the sensitivity to displacements), the now time-dependent photon pressure will induce additional significant oscillations. In contrast to a constant coupling, this effect cannot be canceled by preparing the mechanics in an appropriate initial state. However, by adding a time-dependent linear potential term of the formĤ ext = ω 0 k(τ ) N a (b † +b), all contributions of N a to x m cancel. While addingĤ ext does not modify the QFI for the measurement of displacement, the potential must be known to the same precision as the gravitational field that is being measured. We conclude that, given the Hamiltonian (3), the strongest bound is given by the standard deviation ∆x m .

C. Scaling of the sensitivity given the cavity field restrictions
From the expressions in (41) and (42), we see that increasing the photon number standard deviation ∆N a decreases the spread ∆g 0 . However, since ∆N a must obey the restrictions we derive above, and since these restrictions scale with time, we can consider the fundamental scaling of the QFI when the photon number restriction is taken into account. We focus specifically on the scaling with n, which is a positive integer given by τ = 2πn.
Starting with resonant gravimetry, we identified the requirement that ∆N a l/(2x 0 k 0 ). Since ∆N a does not increase with time, the overall scaling of the sensitivity goes as ∆g 0 ∝ n −1 , where n = τ /(2π), as per (41).
For a modulated optomechanical coupling, we identified the following restriction: ∆N a l/(x 0 k 0 τ sep ). Since ∆g 0 ∝ s −2 ∆N −1 a , as per (42), where s 1, we see that the overall scaling of the sensitivity with respect to s is given by ∆g 0 ∝ s −1 .
These considerations show that a scaling of the sensitivity ∝ τ 2 can be achieved using the modulated coupling, however if the restrictions to the cavity field parameters are taken into account, the scaling is ∝ τ . It remains to be determined whether these restrictions can be circumvented and how they scale with τ when decoherence is taken into account.

D. Limitations due to the Casimir effect
When two objects are placed in close proximity, they will almost always experience a force due to the Casimir effect [88]. While there is an ongoing effort to derive simple expressions for alternative configurations [89], here we use an analytic formula for the acceleration due to the Casimir effect between two homogeneous perfectly conducting spheres, which is given by the spatial derivative of equation (21) of [90] divided by the mass: where c is the speed of light, m is the mass of the sphere, R is the radius and r is the distance between the source and the probe. In our case, the two systems are unlikely to be made of a perfectly conducting material, and they might also not be entirely spherical, but we use (46) to estimate the order-of-magnitude of the resulting Casimir-Polder effect.
In Section VI A, we estimated the fundamental sensitivity of an optomechanical system to the gravitational field produced by a small oscillating sphere. Assuming that both the source mass and the optomechanical probe system are made of tungsten, and that they both weigh ∼ 200 ng, we find an acceleration of the order of ∼ 1 × 10 −12 m s −2 due to the Casimir effect at a distance of 100 µm. The constant gravitational acceleration from the same system is also of order ∼ 1 × 10 −12 m s −2 . This shows that the Casimir effect can become an important systematic factor for gravimetry in the regime that we are considering.
The numbers shown here can be reduced significantly by considering larger distances, or by using a material in-between the source mass and the sensor that acts as a shield to the Casimir effect [91,92]. The addition of the shield induces a stationary Casimir force and the only remaining time-dependent force on the sensor will be the oscillating gravitational field. Here, measuring oscillating gravitational fields instead of static ones has a clear advantage. The only limitation is the size of the shield itself. Additional reductions of the Casimir force can be achieved by adding nanostructures to a metallic surface [93], compensating or modulating the Casimir force with radiation pressure [94] or optical modulation of the charge density [95]. Theoretical investigations also indicate that its sign can be inverted with a shield made out of a left-handed metamaterial [96].
A specific version of the shielding scheme arises in levitated optomechanics when the oscillating source mass can be placed behind the end mirror of the cavity. Then, the mirror itself serves as a shield for Casimir forces [68].

VIII. CONCLUSIONS
In this work, we computed the fundamental sensitivity for time-dependent gravimetry with a nonlinear optomechanical system. We considered both coherent and bright squeezed states of light, and we found that it is possible to significantly enhance the sensitivity of the system by modulating the optomechanical coupling. To ensure that these sensitivities are not influenced by the initial state of the mechanical element and can be achieved through measurements of the cavity state, we identified the points at which the mechanical oscillator and optical mode evolve into a separable state. In addition, we proved that for coherent states the QCRB is saturated for homodyne measurements when the optical mode and mechanical oscillator are found in a separable state. For squeezed coherent states, we found that this is also true when the vacuum contribution is negligible.
Our results serve as a proof-of-principle that an optomechanical system could potentially be used to measure the gravitational field from oscillating source masses as small as 200 nano-grams at a distance of 100 µm. We also provide bounds for quantum optomechanical systems in the nonlinear regime when used as gravitational wave detectors. To successfully detect passing gravitational waves, we have assumed parameters that are experimentally challenging to implement, but not beyond the reach of technological advancement. Finally, it should be noted that our methods can be extended to additional experimental platforms, as long as the Hamiltonian is of the general form (3).
where we have moved into a frame that rotates with the light, and where the operators are defined aŝ Furthermore, the dynamical F coefficients in (B2) are given by where the complex function ξ is given by and where α and β are Bogoliubov coefficients given by Here, P 11 and I P22 are solutions to the following differential equations: P 11 + (1 + 4 D 2 (τ )) P 11 = 0 , with the initial conditions P 11 (0) = 1 andṖ 11 (0) = 0 and I P22 (0) = 0 andİ P22 (0) = 1. Furthermore, the J coefficients in (B2), which arise from the inclusion of a modulated mechanical frequency are given by The full derivation of these quantities and additional examples of their use can be found in [27]. We now present solutions to the F and J coefficients for different choices of k(τ ), D 1 (τ ) and D 2 (τ ). We focus on displaying FN aB± and FB ± , since these make up the QFI in (14) and are used elsewhere both in the main text and in other appendices. While the QFI also depends on FN a , this is often an extremely long and cumbersome expression, and we do not print it here.

Dynamics for a modulated coupling and time-dependent gravitational field
For the case of a modulated optomechanical coupling k(τ ) = k 0 cos(Ω k τ + φ k )) and a time dependent gravitational field D 1 (τ ) = −d 1 (a + cos(Ω d1 τ + φ d1 )), and no modulation of the mechanical frequency, such that D 2 (τ ) = 0, we obtain and the J coefficients in (B8) are again given by J b = τ , J ± = 0.

Dynamics for a time-dependent gravitational field and modulated mechanical frequency
For a constant optomechanical coupling k(τ ) = k 0 , a time-dependent gravitational field D 1 (τ ) = −d 1 (a+ cos(Ω d1 τ + φ d1 )) and a mechanical frequency that is modulated with D 2 (τ ) = d 2 cos(2τ + φ d2 ), we must first find the approximate resonant solutions of the differential equation (B7). For details on how these solutions can be found see Appendix E in Ref [27]. In short, driving the system at Ω d2 = 2 causes the differential equations (B7) to take the form of Mathieu's equation [97]. It has the following form where the solutions y(τ ) correspond to P 11 and I P22 shown in (B7). We now briefly recap the perturbation theory used in [27] to derive solutions for d 2 1. We define a slow time scale X = qτ , as well as the parameter q = −2d 2 . The solutions y can be taken to depend on both scales, such that y(τ, X). The absolute derivative d/dτ in (B11) can then be split into two independent parts: which means that Mathieu's equation (B11) becomes We then expand the solution y(τ, X) for small q as y(τ, X) = y 0 (τ, X) + q y 1 (τ, X) + O(q 2 ) and insert this into the differential equation above. We first recover the regular harmonic oscillation equation for y 0 , which is the limiting case as q → 0: To solve this equation, we propose the following trial solution: Here, A(X) is still undetermined. We continue with the equation for y 1 . To first order in q, we find Inserting our solution for y 0 , we find This expression can be rearranged into We expand the exponentials to find In order for the solution to be stable, we require that secular terms such as resonant terms e iτ vanish. If these do not vanish, the perturbation y 1 will grow exponentially [97]. We also neglect terms that oscillate much faster, such as e 3ix . This leaves us with the condition that which can be differentiated again and solved with the trial solution A(X) = (c 1 − i c 2 ) e (X+iφ d2 )/2 + (c 3 − i c 4 ) e −(X−iφ d1 )/2 for the parameters c 1 , c 2 , c 3 and c 4 . From the requirement in (B20), it is now possible to fix two of the coefficients in (B22). We differentiate A(X) and use (B20) to find that the conditions c 1 = c 2 and c 3 = −c 4 must always be fulfilled. Therefore, A(X) becomes We then recall that X = qx and after combining some exponentials, we obtain the full trial solution for the zeroth order term y 0 : Using the fact that q = −2d 2 , and rearranging, we find The coefficients are then fixed by the initial conditions, which for P 11 read y 0 (0) = 1 andẏ 0 (0) = 0, and I P22 read y 0 (0) = 0 andẏ 0 (0) = 1. Using these, we find the following solutions Using (B4), we can derive the F coefficients. They are however rather lengthy, so we will not display them here.
usual squeezing operator and ζ = re iϕ , though one can readily move between definitions using the standard braiding relations. By (19) the relevant quantity is (∆N a ) 2 |µc,ζ . This can be calculated in a number of ways, but a convenient approach is to first define a new operatorâ ζ as the linear combination [98], where we adopt the usual convention u = cosh(r) and v = e iϕ sinh(r) satisfying |u| 2 − |v| 2 = 1. Then |µ c , ζ are eigenstates ofâ ζ ,â Similarly, we can transformâ ζ back toâ through,â It is also useful to note the commutation relations, With these expressions it is then straightforward to show, where we have used (C12a) to passâ ζ throughâ in the second line. Higher order terms can be found in a similar manner, leading to the following useful results, Using (C13) and (C17) one can then show that the photon number variance is given by Note, in the limit of zero squeezing, we recover |µ c | 2 as expected. A more convenient form is to substitute back in for u and v. With a little algebra we find, Since the photon number variance enters into the QFI (see (C8)), we note that the QFI is maximised when e iϕ 2 µ c is purely imaginary. We then have an enhancement proportional to e 4r over coherent state driving of the mirror (for large photon number), along with a squeezing dependent vacuum contribution.

Resonant driving and modulated coupling
If we assume that k(τ ) = k 0 cos(Ω k τ + φ k ) and D 2 = 0. Again, ξ = e −iτ . If the modulation of the coupling and the oscillations of the gravitational field are at the same frequency, i.e. Ω := Ω d1 = Ω k , the QFI becomes For complete resonance with the mechanics, i.e. for Ω = Ω d = Ω k = 1, the QFI reduces to sech(2r T ) 4(sin(τ )(2a + cos(τ + φ d1 )) + τ cos(φ d1 )) 2 We see that the scaling with time of the QFI depends on the choice of the phases φ d1 and φ k . In the case of φ d1 − φ k = π/2 for example, the QFI contains terms proportional to τ 4 . This is a highly unusual scaling: Normally, under the conditions that coherence is retained, one obtains a scaling of the QFI ∝ τ 2 , as is the case e.g. for a single harmonic oscillator whose frequency one wants to estimate [99], which in itself represents an advantage over the classical scaling ∝ τ . It implies that being able to maintain coherence over long times pays off much more for the optomechanical system with its nonlinear coupling than for a single harmonic oscillator, and suggests to rather reduce the coupling k and increase τ instead in the presence of decoherence, rather than trying to make the coupling as strong as possible, as this is expected to reduce the coherence time due to enhanced non-classicality.
In Appendix D, we show that light and mechanics disentangle at times that are multiples of sπ, with s integer for fractional frequencies Ω = Ω frac = 1 + 2n 1 /s, where n 1 > −s/2 is an integer. At the decoupling times τ = q sπ, with q integer the QFI becomes

Resonant driving and modulated mechanical frequency
We here assume that the optomechanical coupling is constant k(τ ) = k 0 and that the mechanical frequency is modified as D 2 (τ ) = d 2 cos(2τ + φ d2 )). The general QFI is a long expression that we will not give here. Instead, we refer to the main text for plots and simplified expressions for special cases. In particular, for the case of purely oscillating gravitational fields a = 0 and r T 1, the QFI is maximized for φ d2 = −π/2 and φ d1 = 0 and becomes approximately + e d2τ − 1 (9 cos(2τ ) + 3) sinh(d 2 τ ) + 6 sin 2 (τ ) cosh(d 2 τ ) + 16 cos 3 (τ ) − 1 . and the fractional frequencies Ω frac = 1 + 2n 1 /s, where n 1 > −s/2 to obtain positive frequencies. There are infinitely many fractional frequencies Ω frac for which there exist times that are multiples of π at which light and mechanics decouple. Furthermore, from the structure of Ω frac , we see that for each multiple of τ sep given by τ = q sπ, we can find anñ 1 = q n 1 such that 1 + 2ñ 1 /(qs) = Ω frac . Therefore, |KN a | 2 vanishes for all times that are multiples of sπ. For a given Ω frac , the smallest decoupling time sπ is given by the smallest integers n 1 and s > 0 whose quotient n 1 /s is equivalent to (Ω frac − 1)/2. case h(x, ψ τ ) = f 2 (x, ψ τ ) (where we adopt the convention that the hat onx is post squaring, unless the power is after the argument). If the cavity state is initially in a coherent state, this can be achieved when FN 2 a is engineered to be a multiple of 2π (at the decoupling time), then 7 , i.e., |μ c is again a coherent state. In [25] it was shown that for this special case Homodyne measurements saturate the QFI for a constant gravitational field, while in [26] this was confirmed numerically when the parameter of interest is also encoded in a constant frequency shift (along with the observation that Heterodyne measurements preserve a similar scaling). Here we provide an alternate derivation which holds for arbitrary modulations, provided the measurements are performed when the optical and mechanical modes completely disentangle.

Homodyne measurements
We begin by computing the CFI for homodyne measurements.

a. Coherent states
For Homodyne measurements the relevant POVM is constructed via the state |x = |x λ , defined as the eigenstate of the operator,x λ =â A rearrangement of its eigenvalue equation leads to the following action ofâ as, To calculate the R term, we note that under the restriction that |ψ τ = |μ c , then from (E6) and (E9) it is straightforward to show, The last step is to take the expectation value of this function once x λ is promoted to an operator. Here it is useful to note that 8 where H n are the Hermite polynomials. Finally, with a little algebra we find, Note, whenμ c e −iλ is purely imaginary, the CFI simplifies to, which is exactly the QFI given in (19). Note however, that whenμ c e −iλ is purely real the CFI is zero. 7 In the frame rotating with the optical field 8 This identity can be easily derived by noting that for any pair of coherent states |µ 1 and |µ 2 , andD is the usual displacement operator. By evaluating the overlap, and using the generating function of Hermite polynomials, e 2xt−t 2 = ∞ n=0 Hn(x) t n n! , we find the desired result.

b. Squeezed states
For squeezed initial cavity states |µ c , ζ =Ŝ(ζ) |µ c , with squeezing parameter ζ = re iϕ , we can evaluate the CFI using a similar approach to above. The first term in (E5) is independent of the POVM and is given by B 2 N 2 a ζ . This can be found immediately from the expectation value (C17). In order to calculate the corresponding R terms, we must evaluate the overlap x λ |N aÛcŜ (ζ)|µ c at the chosen decoupling time. Again, we will consider the special case where FN 2 a is a multiple of 2π, then by using the identity e −AN 2 aâ e AN 2 a = e A(2Na+1)â and considering the action on a coherent state, we find x λ |N aÛcŜ (ζ)|µ c = x λ |N aŜ (ζ)Û c |µ c = x λ |N a |μ c ,ζ , whereζ = re iφ withφ = ϕ − 2FN a . Following Appendix C 1, we define the operatorâζ =Ŝ(ζ)âŜ † (ζ), with forward and backwards transformations given explicitly by,âζ whereâζ |μ c ,ζ =μ c |μ c ,ζ , and the functionsũ = u = cosh(r) andṽ = e −2iFN a v = e iφ sinh(r) depend explicitly on the parameter ζ. A rearrangement of (E14a) (and the adjoint of (E14b) respectively) gives, The next step is to update (E9) to remove the explicit dependence onâ. Using the adjoint of (E15b) and collectinĝ a † terms on the left hand side, we have, Similarly the overlap, Expanding, and making use of the (tilded) commutation relation (C12b) we find, In order to calculate the term R = h(x λ , ψ τ ) ψτ we need the expectation values ofx λ on the squeezed states |μ c ,ζ (up to fourth power). Here we note that the operatorx λ can be written in terms ofâ ζ andâ † ζ as, where w =ũe iλ −ṽe −iλ ≡ |w|e iλ . Thus we can see that the required expectation values can be found in a similar way as to those over coherent states, μ c ,ζ|x n λ |μ c ,ζ = |w| n (2i) n H n i μ c ,ζ|x λ |μ c ,ζ |w| .
Together with (C17) we can now evaluate the CFI exactly. However, one can greatly simplify the problem by noting that the aim is to find the optimal bound. From (C19) we can expect that the maximum CFI also occurs when R[e − iϕ 2 µ c ] = 0, while in the limit of zero squeezing, (E12) suggests the condition R[μ c e −iλ ] = 0. Note, in the latter one needs to choose λ based not only on the initial state, but on the additional phase FN a picked up after the evolution. Thus the optimal λ varies with time. Writing µ c = |µ c |e iχ these two conditions imply, and so w = e −r e iλ , with x λ = 0. With these simplifications, and together with (C17) for the first term in (E5), the CFI is given by, This is less than the maximum QFI by only a vacuum contribution. In general, however, the CFI can still be non-zero when µ c = 0, I (hom) ζ (µ c = 0) = B 2 2 sinh 2 (2r) sin 2 (φ − 2λ) (cosh(2r) − sinh(2r) cos(φ − 2λ)) 2 . (E22) This reaches a maximum of I hom ζ = 2B 2 sinh 2 (2r) forφ − 2λ = ±2[nπ ± tan −1 (e −2r )]. However, for all but very small photon number (and large r) the optimal CFI is given by (E21).

Heterodyne measurements
The Heterodyne measurement case is somewhat more straightforward as the probabilities are calculated with respect to coherent states. Note that an additional factor of 1/π appears in the definition of I, (E1), to account for the identity operator in the coherent state basis though this is then removed when moving to the expectation value expressions.

a. Coherent states
When there is no squeezing of the initial cavity states the CFI can be evaluated quickly. Replacing |x = |β in (E6) we have, and so, The CFI is then given by, which is exactly half of the QFI (19). Thus the precision obtainable through Heterodyne measurements also preserves the optimal scaling behaviour, though in general is less favourable than performing Homodyne measurements.
Note, e − iφ 2μ c = e − iϕ 2 µ c , which means one can fix the third term to zero (and thereby maximise I het ζ ) by choice of the initial state alone. However, even for large photon numbers, the CFI is smaller than the QFI by a factor of 2e r cosh(r).

Mechanical expectation values
In this section, we compute x m for the different cases of the dynamics considered in the main text.
The conditions that follow from the above analysis are that 2 x 0 k 0 N a l and 2 x 0 k 0 ∆N a l such that the interaction Hamiltonian is still valid.
The conditions that follow from the above analysis are that 2 x 0 (k 0 N a +d 1 a) l, x 0 d 1 τ l and 2 x 0 k 0 ∆N a l such that the interaction Hamiltonian is still valid. We conclude that the restrictions on N a and ∆N a do not increase with τ . Furthermore, the driving D 1 does not affect the restriction on the standard deviation ∆N a as it does not change its evolution.
Finally, we plot the expectation value x m and the standard deviation ∆x m as a function of time τ in Figure 5. quantities for resonant gravitational fields, a jointly resonantly modulated coupling and gravitational field, jointly modulated coupling and gravitational field at the fractional frequencies identified in Appendix D, and jointly modulated mechanical frequency and gravitational field, respectively.