Gravimetry through non-linear optomechanics

Precision gravimetry is key to a number of scientific and industrial applications, including climate change research, space exploration, geological surveys and fundamental investigations into the nature of gravity. A variety of quantum systems, such as atom interferometry and on-chip-Bose–Einstein condensates have thus far been investigated to this aim. Here, we propose a new method which involves using a quantum optomechanical system for measurements of gravitational acceleration. As a proof-of-concept, we investigate the fundamental sensitivity for gravitational accelerometry of a cavity optomechanical system with a trilinear radiation pressure light-matter interaction. The phase of the optical output encodes the gravitational acceleration g and is the only component which needs to be measured. We prove analytically that homodyne detection is the optimal readout method and we predict an ideal fundamental sensitivity of Δg = 10−15 ms−2 for state-of-the-art parameters of optomechanical systems, showing that they could, in principle, surpass the best atomic interferometers even for low optical intensities. Further, we show that the scheme is strikingly robust to the initial thermal state of the oscillator.

T he practise of measuring the gravitational acceleration galso known as gravimetry-has led to important advances in both fundamental science and industry. For example, local gravity variations due to mass redistribution driven by climate change have been mapped with the GRACE satellite [1][2][3] , and more recently, the Juno spacecraft mission reported the measurement of the gravity harmonics of Jupiter 4 . Furthermore, precise measurements of g can test for small deviations from Newtonian gravity on extremely small scales, which may provide indications of a deeper theory of quantum gravity 5 . In industry, precision accelerometry is extensively used in inertial navigation technologies and for conducting geological surveys.
While classical systems have long been utilised to perform accurate measurements of g, quantum systems offer several useful advantages, including reduced noise levels, a compact setup and most importantly an increased measurement sensitivity achieved through the power of coherence and interferometry. Over the past decade, a variety of quantum systems have been explored to this aim, in both theory and practice. The largest research effort to date has focused on atom interferometry [6][7][8][9] , for which the highest achieved sensitivity currently stands at Δg = 4.3 × 10 −9 ms −29 . A similar investigation has been carried out for both onchip and fountain Bose-Einstein condensate (BEC) interferometry with best sensitivity Δg = 7.8 × 10 −10 ms −210 . Finally, a proposal for using magnetically levitated spheres that predicts sensitivities of 2.2 × 10 −9 ms −2 Hz −1/2 has been put forward in 11 . For comparison, the current commercial standard is set by the LaCoste FG5-X gravimeter which can achieve a measurement sensitivity of 1.5 × 10 −9 ms −2 Hz −1/212 . More generally, the broader topic of using quantum systems to probe relativistic phenomena is currently being pursued with great interest (see for example [13][14][15][16][17][18][19][20].
A key advantage to quantum systems are their interferometric properties. The following question arises: How can these interferometric properties be enhanced in order to improve the measurement sensitivity? One possibility is to place a quantum system in the form of a mechanical oscillator in an optical cavity, a research area known as quantum optomechanics 21 . See Fig. 1 for an illustration of a nanodiamond trapped in an optical cavity as an example of a class of optomechanical systems. The addition of the cavity allows for a strong coherent coupling between light and oscillator which, as we shall see, cancels out any initial thermal noise and fundamentally improves the measurement sensitivity of the device.
Within classical optomechanics, the idea of gravimetry and accelerometry by optically detecting the mechanical oscillator has been experimentally realised by Cervantes et al. 22 . Other avenues, such as the detection of high-frequency gravitational waves through the driving of resonant mechanical elements was proposed also in 23 . In the related field of electromechanics, Schrödinger cat states and a Kerr nonlinearity have recently been found to be useful for the same applications 24 . However, the ensuing fundamental limits on the measurement sensitivity of gravimetry in the quantum regime of optomechanics using its trilinear radiation pressure interaction is yet to be investigated. Here we undertake this task and obtain some striking results: Firstly, it is possible, in principle, to surpass the sensitivity Δg that has been obtained in atom interferometers and other implementations to date. Secondly, due to the periodic decoupling of light and mechanics, the mechanical element does not require initial cooling to the ground state to improve the fundamental sensitivity of the gravimeter and, finally, the best possible sensitivity is achieved by a simple homodyne measurement of the cavity field, while only a low photon number in the cavity is required. That is, no measurement on the mechanical oscillator is required. Unlike the case of atomic interferometers, in optomechanics the interaction of light and matter is continuous, and we will see that our Hamiltonian cyclically entangles and disentangles the light and mechanics, leading to their decoupling. It follows that the experimental challenge will be to maintain the quantum coherence of the field and mechanics over the duration of each run of the experiment, which we set as one oscillation period of the mechanical element. This requirement, on which the plausibility of the scheme hinges, will be discussed in some detail.

Results
Optomechanics and Newtonian gravity. Let us begin by considering a general optomechanical system consisting of a mechanical oscillator coupled to a light-field in the cavity. The non-gravitational Hamiltonian that describes the dynamics of an optomechanical system is given by: 25,26 H ¼ hω câ yâ þ hω mb yb À hkâ yâ ðb y þbÞ; ð1Þ whereâ;â y are the annihilation and creation operators for the cavity field with frequency ω c ,b;b y are the annihilation and creation operators for the mechanical oscillator with frequency ω m , and k (usually denoted g in the literature, but this we shall here reserve for gravity) is a coupling constant that determines the interaction strength between the photon numberâ yâ and the positionx o / ðb y þbÞ of the oscillator. In order to introduce a coupling to a gravitational potential in the Hamiltonian, we add a term of the form mgx o cos θ. Here, m is the mass of the mechanical oscillator, g is the gravitational acceleration,x o ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h=2mω m p ðb y þbÞ is the position operator acting on the mechanical oscillator and θ is an angle from the vertical axis that we include in order to describe inclined systems, similar to 27 . Note that while the mass m appears as a coupling in the Hamiltonian, we will later see that measurements of g are mass-independent, which is what we expect from the equivalence of inertial and gravitational mass. With the addition of Newtonian gravity, the Hamiltonian of the system thus becomeŝ where L is the length of the cavity and m is the mass of the mirror. A levitated nano-or micro-crystal (e.g. a diamond or silicon bead), on the other hand, has a k given by 28,29 where ε 0 is the permittivity of free space, V c is the cavity mode volume, and k c is the wavevector of the laser, given by 2π/λ, where λ is the laser wavelength. P = 3Vϵ 0 (ϵ -1)/(ϵ + 2) is the polarizability of the levitated object of volume V and ϵ is the relative electric permittivity. Alternatively, we can also consider a BEC trapped in a cavity. Here, the collective motion of the ensemble acts as the massive oscillator. For this system, the coupling constant is given by 30,31 where N is the number of atoms in the ensemble, g 0 is the singleatom cavity QED coupling rate, M = Nm is the collective mass of all the trapped atoms with individual mass m, k l is the wavevector of the laser and Δ ca = ω p -ω c with pumping frequency ω p . We will return to these expressions when computing the fundamental sensitivity limits for each system in the latter part of the paper.
System dynamics. In order to simplify the time evolution operatorÛðtÞ corresponding to the above Hamiltonian, we rescaleĤ G by dividing all terms by the oscillator frequency ω m . As a result, the time parameter t now represents the labframe time multiplied by ω m , such that the oscillator has undergone a full oscillation cycle at t = 2π. The operatorÛðtÞ can then be written in the following decoupled form (see ref. 26 for details of the derivation in the absence of gravity): where . As a rule, we will denote any dimensionless quantity with a bar. For time-dependent variables, such as dissipation rates, this means they have been rescaled with respect to ω m .
We now assume that the cavity field mode and the mechanical oscillator are initially in coherent states |α〉 C and |β〉 O respectively. For laser light injected into the cavity, this is the natural assumption. The oscillator, on the other hand, will in reality be initialised as a thermal state, which corresponds to a random coherent state |β〉 O according to a thermal distribution. However, by starting out with a coherent state we will later argue that the gravimetric phase accumulated by the light does not depend on jβi O so that our procedure works equally well for an arbitrary thermal state. A formal proof of this statement can be found in Supplementary Note 1. The initial state at t = 0 is then given by jΨ 0 ð Þi ¼ jαi C jβi O , and underÛðtÞ it gives us the following state where τ = tsin t, and |ϕ n (t)〉 O are coherent states of the oscillator given by ϕ n t ð Þ O ¼ e Àit β þ kn À g À Á 1 À e Àit ð Þ . In the derivation of this state, we have adopted a rotating frame for the cavity field, thus ignoring the free evolution induced by the term exp{-ira † a}.
The state in Eq. 7 show us that light and mechanics will entangle and disentangle periodically, with maximum entanglement occurring at t = π. At t = 2π, the oscillator state ϕ n O returns to β j i O regardless of the values of k; g and β, and therefore by extension a thermal state also returns to its initial state because it will undergo the same compact evolution. This means that the initial oscillator state does not impact the fundamental sensitivity of this scheme. As already mentioned, a formal proof of this can be found in Supplementary Note 1. Most importantly however, at t = 2π the cavity state is completely decoupled from the oscillator, meaning that all information about g is transferred to the phase of the cavity state. As a result, any measurement scheme only needs to consider the cavity state after one oscillation period, meaning that direct or indirect access to the oscillator state is not required. This will greatly simplifies an experimental implementation, as measuring the oscillator state is generally difficult. This convenient property arises from the interferometric properties of the oscillator; its quantum nature allows it to acts as an interferometer to ensure that any initial thermal noise is removed from the cavity field, and thereby our scheme does not require cooling of the oscillator to a pure ground state. In other words, our results are valid for both coherent and thermal states. Note however that decoherence ensuing from damping to the oscillator motion during the state evolution will adversely affect the final measurement sensitivity and cause the oscillator state to grow increasingly mixed. We will not consider this kind of decoherence in this work, and instead assume that the mechanical element remains coherent over one oscillation period.
We can visualise some of the dynamics of the state in Eq. 7 by computing the expectation values of the field quadraturesX c ¼ ðâ y þâÞ= ffiffi ffi 2 p andP c ¼ iðâ y ÀâÞ= ffiffi ffi 2 p 32 . We focus on the cavity state, which we obtain by tracing out the oscillator. The tracedout cavity state is given by For decoherence-free evolution, the trajectories traced out by the system in phase space can be seen for different values of g in Fig. 2a with k ¼ g ¼ 1 and Fig. 2b with k ¼ 1; g ¼ 2. Both figures are plotted with α = β = 1. We observe that the system performs increasingly complex trajectories for larger values of g.
We noted above that the light and mechanics periodically entangle and disentangle during its evolution. In order to see this more clearly, we can compute the linear entropy S(t) for the traced-out cavity state ρ c (t) in Eq. 8. The linear entropy is defined as The linear entropy tells us about the entanglement between the cavity and oscillator states. The results can be found in Fig. 3a for pure state and in Fig. 3b for states undergoing decoherence with photon dissipation rate κ ¼ κ=ω m . We see that S(t) increases until the state is maximally entangled at t = π. While a pure state completely decouples the light and mechanics at t = 2π for any values of k and g, a decohering state becomes increasingly mixed and does not return to its original state.
Quantum metrology. We now come to our main results which concern the use of optomechanical systems as gravimeters. The question we wish to answer is: what is the best fundamental sensitivity Δg with which an optomechanical system can measure the gravitational acceleration g? Here, Δg denotes the standard deviation of a gravimetric measurement. We can directly predict Δg from the system's dynamics by calculating the Fisher information I F (t) which provides a natural lower bound on the variance Var(g) of an unknown parameter, in our case g. This relationship is captured by the Cramér-Rao inequality [33][34][35] Var where N is the number of measurements. Thus if we maximise I F (t), we minimise the measurement spread of g.
Quantum Fisher information. The Fisher information comes in two forms: the measurement-specific classical Fisher information (CFI) and the quantum Fisher information (QFI). The QFI, which we denote H Q (t), is computed by optimising over all possible positive-operator valued measures (POVMs) and their resulting CFI 36 . Thus H Q represents the ultimate bound on obtainable information from a system, but it does not reveal which specific measurement is required to achieve it. For a general mixed quantum state ρ(t) the QFI is given by where L is the symmetric logarithmic derivative. In general, it is difficult to obtain L analytically, especially for noisy systems. There are however methods for finding a noisy bound on the Cramér-Rao inequality 37 . A similar method for many-body systems was proposed in 38 , and numerical methods were shown to be effective for a class of specific systems 39 . We shall not be using these methods here, as we shall instead investigate specific measurements for the noisy scenario to better  approximate an experimental setting. This will later allow us to prove the optimality of the homodyne measurement. Let us start by deriving a fundamental bound to the sensitivity. We specialise to the simpler case where the state ρ(t) is pure. Setting ρ(t) = |Ψ(t)〉〈Ψ(t)| the QFI becomes where we have used the notation ∂ g = ∂/∂ g . At first glance, the QFI of the global system might not seem very relevant as the mechanical part of the optomechanical system cannot easily be measured directly. However, we recall that the coherent state ϕ n t ð Þ O returns to β j i O at t = 2π, so that all information about g is transferred to the phase of the pure, decoupled cavity state. Since the decoupling time does not depend on β, this is also the case for a thermal state that may be written as a statistical mixture of coherent states (see Supplementary Note 1 for a proof of this statement). Calculating the QFI for this state will therefore provide an experimentally accessible notion of the fundamental sensitivity of the device. We find the following expression for H Q (t) at t = 2π: Note that the mass term m is cancelled by the appearance of √m in the coupling constant k, so that the final accelerometry measurement will be mass-independent. We also note the strong dependence on k and ω m , and that the expression scales linearly with the number of photons α j j 2 . To find H Q (t) for the global state at any time t we resort to numerical calculations. We consider the case k ¼ g ¼ 1 to allow for future comparisons with subsequent numerical evaluation of the CFI, which will be restricted to the same narrow parameter range. The resulting H Q (t) as a function of t can be found in Fig. 4a with k ¼ g ¼ 1 and β = 1 for various values of jαj 2 . We note that H Q (t) reaches its maximum value at t = 2π, which means that Eq. 12 returns the largest possible value during one oscillation period for any choice of system.

Classical Fisher information.
Let us now consider a specific measurement of g. The CFI I F (t) determines the minimum standard deviation of a parameter estimator once we have chosen a single specific measurement with POVM elements fΠ x g. The CFI is given by the expression where p xjg ð Þ ¼ tr Π x ρðgÞ ½ is a conditional probability distribution.
We now consider a general homodyne measurement on the traced-out cavity state ρ C . For notational convenience, we use a general Hermitian operatorx λ ¼ ðâ exp Àiλ f gþâ y exp iλ f gÞ= ffiffi ffi 2 p , where λ denotes a label that rotates between the field quadratures 40 . Any two operators that differ by λ = π/2 form a Fisher information for measurements of gravitational acceleration g. Plots of the quantum Fisher information (QFI) and classical Fisher information (CFI) for measurements of g. a shows the QFI H Q (t) for 1, 4 and 9 photons jαj 2 with rescaled couplings k ¼ g ¼ 1. b shows the dimensionless CFI I F ðtÞ with and without photon dissipation rate κ ¼ κ=ω m for a momentum measurement (λ = 1/2) and k ¼ g ¼ 1. c shows how the peak of the CFI for a momentum measurement narrows with increasing k and with constant g ¼ 1. d shows the I F (t) for a momentum measurement of photons that leak from the cavity with environmental coupling γ ¼ 0:1 and k ¼ g ¼ conjugate pair which satisfies the position-momentum commutator relation. In the following, we shall refer to the choices λ = 0 and λ = π/2 as a position and momentum measurement respectively. In order to calculate I F (t) we must find the probability distribution pðx λ jgÞ ¼ tr½jx λ ihx λ jρ c ðgÞ, where jx λ ihx λ j are the eigenstate ofx λ . While the position eigenstates themselves are not proper vectors, we can make use of a standard result from the quantum harmonic oscillator: where H n (x) are the Hermite polynomials of order n. These probabilities in turn gives rise to a CFI of the form ðnÀn′Þc n;n′ d n;n′ ðx λ Þf n;n′ ! 2 P n;n′ c n;n′ d n;n′ ðx λ Þf n;n′ ; ð15Þ where c n;n′ ¼ ðα Ã Þ n′ α n ffiffiffiffiffiffiffiffiffi ffi n!n′! p e ið k 2 ðn 2 Àn ′2 ÞÀ2 k gðnÀn′ÞÞτ ; ð16AÞ Timescales of measurements. Let us analyse the expression for I F (t). We immediately note that any terms in the sum with n = n′ do not contribute to the Fisher information. The remaining behaviour of I F can be inferred from the second exponential in f n, n′ , namely exp Àjϕ n′ j 2 =2 À jϕ n j 2 =2 þ ϕ Ã n′ ϕ n È É as this will dominate the entire expression for large k. If we simplify the expression in the exponential, we find that it is equal to exp À k 2 ðn À n′Þ 2 ð1 À cos tÞ þ kðn À n′Þ For n ≠ n′ and large k, the first term will dominate, and the exponential will be small for any t that is not a multiple of 2π. In other words, the Fisher information for a homodyne measurement becomes significant only when light and mechanics are completely decoupled. Figure 4c shows how the CFI for a momentum measurement (with λ = π/2) for g ¼ α ¼ β ¼ 1 and k ¼ 1; 2; 5 becomes increasingly narrow as k grows larger. For clarity, we have rescaled I F with k in the plot. Note that for small k we still find large I F at times t ≠ 2π.
We saw earlier that the QFI scales with k 2 , which mean that the scheme favours systems with a large single-photon coupling. We shall soon show that the CFI coincides with the QFI at t = 2π, but in the meantime we must explore what the narrowing of the CFI at t = 2π entails for our measurement scheme. The narrow peak of the CFI will require the homodyne measurement to be performed within an increasingly narrow time-window. We can estimate the timescale in question by finding the full-width-halfmaximum (FWHM) of the peak. To do so, we consider only the dominant first term À k 2 ðn À n′Þ 2 ð1 À cos tÞ for small perturbations in t around t = 2π, thus cosð2π þ t ′ Þ % 1 À t ′2 =2. That brings the first term into the form À k 2 ðn À n′Þ 2 t ′2 =2, which is now a Gaussian distribution. For a Gaussian function with exp Àðx À x 0 Þ 2 =ð2σ 2 Þ È É , the FWHM is given by 2 ffiffiffiffiffiffiffiffiffiffi ffi 2 ln 2 p σ. In our case, we find σ 2 ¼ ð2 k 2 ðn À n′Þ 2 Þ À1 . We already noted that terms with (n − n′) will not contribute to the CFI, and any term with jn À n′ ) 1 will just cause the peak to narrow further. Thus we only consider the terms with jn À n′j ¼ 1, leaving us with σ ¼ ð2 kÞ À1 , and so we conclude that any measurement must be performed roughly on a timescale of ðω m kÞ À1 ¼ k À1 .
Optimality of homodyne detection. Let us see if we can simplify the expression for I F (t) even further and whether it bears any semblance to the QFI. At t = 2π, ϕ n (2π) = β and η = 0. Then setting k and g to integer values causes I F (2π) to lose all dependence of g. The coefficients reduce to c n;n′ ¼ ðα Ã Þ n′ α n = ffiffiffiffiffiffiffiffiffi ffi n!n′! p and f n;n′ ¼ 1. We now consider the generating function for the Hermite polynomials e 2xtÀt 2 ¼ P 1 n¼0 t n H n =n!. Taking the derivative of both sides results in ð2x À 2tÞe 2xtÀt 2 ¼ P 1 n¼n′ t nÀn′ H n =ðn À n′Þ!, which we can use to show that Eq. 15 reduces to the compact expression This expression coincides precisely with the QFI in Eq. 12 for complementary choices of λ and α. To better see why, we rewrite the term in the brackets as ðe Àiλ À e iλ Þi Re α f g À ðe Àiλ þ e iλ ÞIm α f g Â Ã 2 .
We now note that when λ = 0, only Im α f g contributes to the CFI, whereas at λ = π/2, only Refαg contributes. For both of these specific choices of λ, and when matched by α being either entirely real or imaginary, the CFI coincides precisely with the QFI in Eq. 12 because the term in the brackets reduces to 4Refαg 2 or 4Im α f g 2 , respectively. We conclude that the homodyne measurement saturates the QFI limit up to a phase dependence of α, which can always be accounted for by changing the quadrature of the homodyne measurement. Note, however, that at other times than t = 2π. the homodyne measurement will be zero for all choices of λ and α, and so it only saturates the QFI when the light and mechanics have decoupled.
Finally, the absence of g from I F (2π) is not a problem for sensing g-it just means that the sensitivity at times t = 2π is independent of the actual value of g. Numerical analysis suggests that larger values for g causes the CFI to oscillate increasingly quickly before reaching its maximum value (see Supplementary Note 2). The optimality of the homodyne detection for sensing within our scheme is greatly advantageous as it is a routine measurement which is easy to accomplish. It has in fact also been shown to be an optimal measurement 41 in other contexts.
Decoherence. The calculation above is valid for pure states, but in practice every measurement will suffer various forms of decoherence. We will here investigate the effects of decoherence on the CFI for a narrow parameter range, as realistic parameters are very difficult to simulate numerically. We shall later use these results as indications of the behaviour of realistic systems.
There exists a large variety of decoherence effects for optomechanical systems, such as decoherence due to photons leaking from the cavity, or phonons gradually being lost from the mechanical element. The latter manifests as a gradual damping of the oscillator motion, which moves the state towards a mixture in the coherent state basis [42][43][44][45] . This problem has previously been treated analytically, which is possible because the decoherence operators commute with the Hamiltonian. Thus we refer to these works and will not treat the mechanical decoherence here. Instead, we make the assumption that the phonon decoherence is negligible over one oscillation period of the oscillator.
The effect of photons leaking from a cavity on a state ρ(t) can be modelled using a Lindblad master equation of the form ∂ρðtÞ ∂t where fÁ; Ág denotes the anti-commutator,L ¼ ffiffi ffi κ pâ are Lindblad operators and κ ¼ κ=ω m is the decoherence rate κ with respect to the rescaled time t. This equation cannot easily be solved analytically since the operatorâ does not commute with the HamiltonianĤ G in Eq. 2. Some solutions have been found for specific cases, for example when assuming that the photon leakage occurs only during the injection of the state into the cavity. The decoherence can then be modelled as a series of beamsplitters 46 . We will not consider these modifications here, but instead solve the Lindblad master equation numerically and compute the Fisher information I F (t) for the resulting mixed state.
In all subsequent numerical evaluations, we will set k ¼ g ¼ 1 and α = 1 (note the choice of α 2 R, which will optimise the CFI for λ = π/2). Larger values will cause the system to quickly grow numerically unstable due to the inclusion of non-linear terms such as (a † a) 2 in the evolution in Eq. 6. While k ¼ 1 is experimentally achievable with the right choice of parameters, we can justify setting g ¼ 1 by noting that it physically corresponds to a heavily inclined cavity with θ ≈ π/2. Since we are interested in the general behaviour of the CFI under decoherence, we will here be working with the dimensionless Fisher information I F ðtÞ (see the Methods section). Thus these numerical investigations should only be seen as a indication as to how decoherence will affect I F (t), and not as predictions for the sensitivity of a realised device. We shall later extrapolate from these results to make a prediction for realistic systems.
The behaviour of I F ðtÞ can be found in Fig. 4b for a momentum measurements with λ = π/2. A measurement with λ = 0 and α = i would show the same results. We note that larger values of κ do affect the CFI adversely, but setting κ ¼ 0:1 implies that about 10% of the pure-state CFI is still accessible.
Measurements of leaking photons. In practise, a homodyne measurement is performed by monitoring and measuring the photons that continuously leak from the cavity. Aside from the experimental considerations, such a scheme also negates part of the photon dissipation considered above. We briefly estimate the CFI obtained through such a setup by using a simplified model where a pure vacuum state of the environment |0〉 E is added to our original state |Ψ(t)〉 CO , giving us the combined initial state jΨðtÞi CO j0i E . We then add a rotating wave interaction termĤ I to the HamiltonianĤ G in Eq. 2, of the form where γ is the interaction strength andĉ andĉ y are the creation and annihilation operators of the environment. The effect of this interaction Hamiltonian is to couple the cavity state to the environment which causes information about g to slowly leak out from the cavity into j0i E . As before, we evolve the full state for a single-photon jαj 2 ¼ 1 and with parameters k ¼ g ¼ 1. To maximise the CFI, we choose α 2 R and λ = π/2. The results can be found in Fig. 4d for a rescaled coupling strength γ ¼ 0:1, where γ ¼ γ=ω m . As evident from Fig. 4d, we suffer a 10 -2 reduction in the information that can be extracted from the system. Note also that the behaviour of I F (t) for this scenario will most likely also resemble a delta function centred around t = 2π for realistic parameters. Additional plots for this simplified model can be found in Supplementary Note 3.
Ideal sensitivities. In this section we shall first calculate the ideal Fisher information for the three optomechanical systems considered above, and then discuss the experimental challenges and advantages to an optomechanical gravimeter. As we here calculate the fundamental sensitivity, which is unlikely to be realised, we will only concern ourselves with order-of-magnitude estimates. These results are meant to showcase the potential of optomechanical systems, and to do so we have chosen state-of-the-art parameters that have been implemented in a variety of systems. For discussions of an experimental implementation including noise, see the Discussion.
Starting with the Fabry-Perot cavity system, we choose a fully vertical cavity with θ = 0 and use the following state-of-the-art experimental parameters: We choose a mass m = 10 −6 kg, oscillator frequency ω m = 10 3 Hz, cavity frequency ω c = 10 14 Hz, cavity length L = 10 −5 m and a photon number of jαj 2 ¼ 10 6 . For these values, the rescaled coupling constant in Eq. 3 becomes k FP % 2:30, which gives us a Fisher information of I F = 1.58 × 10 28 m −2 s 4 . This implies a sensitivity of Δg ≈ 7.96 × 10 −15 ms −2 .
Finally, let us also consider cold atoms trapped in a cavity. Based on 30 , we choose the following parameters: a wavelength λ = 780 nm, implying ω c = 10 15 Hz, a single-atom coupling of g 0 = 10 7 Hz, an atomic oscillation frequency ω m = 10 2 Hz, a singleatom mass m = 10 −25 kg, a detuning of Δ ca = 10 11 Hz, and a laser wavevector of k l = 10 8 m −1 . With N = 10 5 atoms trapped in the cavity, we find that k BEC ¼ 2:30 10 6 and I F = 1.58 × 10 19 m −2 s 4 , giving a sensitivity of Δg ≈ 2.5 × 10 −10 ms −2 . The reason for this disparity seems to be that the polarisability of the collection of cold atoms is not high enough to match the polarisability exhibited by the nanosphere. The number of trapped atoms can hardly match the number of atoms in a single nanosphere. One would either have to increase the number of atoms trapped in the cavity or increase the single-atom coupling strength to increase the Fisher information.
Comparison of theoretical results. Let us briefly compare the results obtained here with the performance of other quantum systems in the literature. In Table 1 we have listed a variety of experimentally implemented gravimeter systems with their best achieved sensitivity to date. Table 2, on the other hand, lists the ideal fundamental limits to sensitivities calculated in this work and others. The values for Δg and Δg= ffiffiffiffiffiffi Hz p are presented in units of ms −2 and ms −2 Hz −1/2 , respectively. The last column in Table 1 lists the integration time for each experiment, whereas in Table 2 the last column lists the experimental cycle time set by the oscillation frequency of the system in question. For atom interferometry, it is suggested in 49 that sensitivities of Δg~10 −12 ms −2 might be achieved, and a study of the fundamental limits has very recently been presented in 50 .

Discussion
In this work, we investigated a new scheme for measurements of the gravitational acceleration g using a compact cavity optomechanical system with the usual trilinear optomechanical coupling to the cavity field. We derived a fundamental limit to the sensitivity Δg by computing the QFI and showed that the optimal sensitivity is achieved by a homodyne detection scheme performed on the cavity state at time t = 2π. That is, no direct measurement of the mechanical oscillator is required. Using the expression in Eq. 13 and state-of-the-art experimental parameters, we predict a upper bound on the sensitivity of order Δg1 0 −15 ms −2 for both a Fabry-Perot cavity and a levitated microsphere cavity, and Δg~10 −10 ms −2 for trapped cold atoms. These values compare favourably to all other currently available experimental and theoretical gravimetry proposals (see Tables 1  and 2). Furthermore, the quantum nature of the oscillator ensures that any thermal distribution in its initial state does not affect the fundamental sensitivity. However, as our scheme relies on superpositions involving distinct coherent states, we require thermal decoherence during one period of the oscillator motion to be negligible, which we estimate requires a Q-factor of at least 10 6 for the case of a Fabry-Perot cavity (see below). To explore the effects of photons leaking from the cavity, we numerically explored a narrow parameter range with k ¼ g ¼ 1, which physically corresponds to a nearly horizontally aligned cavity. We found that this form of decoherence does affect the system's performance, but not severely. Finally, we briefly investigated what proportion of Δg we retain by performing measurements on the photons that leak from the cavity. Using a simplified noise model, we found a reduction of 10 −2 in the resulting Fisher information. Given these results, we believe that there is significant potential in the use of quantum optomechanical systems for measurements of gravity and acceleration.
Let us now address some of the experimental challenges related to this scheme. Due to measurement inefficiencies and additional sources of decoherence not considered here, the final performance of optomechanical systems will naturally be expected to be lower than the values presented in Table 2. While we have shown that the initial optomechanical state does not need cooling to the ground state, thermal noise due to external influences during the evolution will gradually decohere the oscillator motion. We estimate that in the case of a Fabry-Perot cavity cooled to a temperature of milliKelvin, a number of ħω m /(k B T th ) = N phonons are present in the system at any time. Here, k B is Boltzmann's constant and T th is the system's temperature. To retain coherence throughout the evolution, we require that κ m N ( ω m , where κ m is the phonon dissipation rate. In other words, the timescale of phonon decoherence κ m must be much less than the characteristic timescale of the system. With ω m = 1 kHz, as we assumed for Fabry-Perot cavities, we find N = 10 5 and κ m = 10 -2 Hz. A cavity which achieves such a decoherence rate must have a mechanical Q-factor of at least Q = ω m /κ m~1 0 6 to retain coherence, a regime which is not unprecedented.
Next, let us discuss which parameters yield the best sensitivities. Firstly, we note that the QFI in Eq. 12 ultimately scales with ω À6 m . In addition to the factor ω À3 m in the denominator, we acquire an extra ω À2 m from the rescaled coupling constant k ¼ k=ω m . The final factor of ω m comes from the dependence of ω m in k 2 . Given this scaling, we require ω m to be as small as possible. At the same time, we also require the photon dissipation rate κ to be low. From our simulations, we saw that we require at least κ ¼ κ=ω m ¼ 0:1. This combination is difficult to achieve as low ω m means the cavity must remain coherent over longer timescales. Therefore, the main experimental challenge of this scheme is to reduce ω m and κ at the same time. Taking our numerical results as guidance, we essentially require that κ ¼ κ=ω m ( 1, which is nothing but the resolved side-band regime 28 . In the above, we used state-of-the-art parameters to calculate the ideal QFI for a variety of systems. However, as we just saw, the photon dissipation rate κ must be very low for these sensitivities to be achieved, and this has not yet been experimentally demonstrated for the parameters we used in Table 2. As technology improves we expect that this to be possible in future experiments, but for now, let us estimate the sensitivities that could be achieved today already. One of the best coherence times to date was demonstrated in 51 , which achieved a cavity linewidth of κ = 660 Hz. To achieve a rescaled photon rate of κ ¼ 0:1 for this system, we let ω m = 6600 Hz and use L = 9.4 cm as reported in the paper. We keep m = 10 −6 kg (since the QFI is ultimately independent of mass) and let ω c = 10 14 Hz as before. Because the oscillation frequency ω m is rather high, we choose to calculate I F for the Fabry-Perot cavity with a mechanical mirror, as this system performed slightly better for higher ω m . The resulting coupling constant is k FP ¼ 1:44 10 À5 , and the Fisher information is I F ≈ 2.16 × 10 15 m −2 s 4 . This leads to Δg ≈ 2.15 × 10 −8 ms −2 . If we now assume that decoherence causes a similar proportion of the Fisher information to dissipate at these parameters compared to the ones chosen in our numerical simulations, we see that we retain about 10% of the pure-state Fisher information. Using this    12 1 × 10 −9 1.5 × 10 −7 6.25 h Atom intf 9 . 5 × 1 0 −9 4.2 × 10 −8 100 s On-chip BEC 10 7.8 × 10 −10 5.3 × 10 −9 100 s Optomech. accel 22 .
3.10 × 10 −5 9.81 × 10 −7 10 −3 s a These include the commercial LaCoste FG5-X, atom interferometry, gravimetry through on-chip Bose-Einstein condensate (BEC) and classical optomechanical accelerometry. The second column lists the sensitivity Δg in ms −2 and the third column lists the ffiffiffiffiffiffi ffi Hz p -noise Δg= ffiffiffiffiffi ffi Hz p in ms −2 Hz −1/2 . The last column indicates the integration time needed to achieve each sensitivity a This value was provided to us by the authors of ref. 22 assumption, we find Δg ≈ 6.80 × 10 −8 ms −2 and a ffiffiffiffiffiffi Hz p -noise of 8.37 × 10 −10 ms −2 / ffiffiffiffiffiffi Hz p . This is directly comparable with the values in Table 1, and so we believe that this scheme could be experimentally realised today, although the experimental challenges are of course substantial.
Let us briefly discuss ways in which we can decrease κ further and how this might affect the Fisher information. A heuristic estimate for κ can be given by considering the number of times per second that a single photon traverses the cavity. Each time the photon is reflected at the mirror, it has a T = 1-R chance of being transmitted instead of reflected. Here, T is the proportion of transmissions and R is the proportion of the number of reflections. The photon bounces off a mirror c/L times per second, where c is the speed of light. Thus we can take the dissipation rate to be κ = Tc/L, which means that increasing L decrease the photon dissipation rate κ, as the photon is effectively spending longer inside the cavity. However, increasing L also decreases the single-photon coupling constant, as we saw in the calculation above. This is true for all couplings we quote here, but it is perhaps most clearly seen for the case of the mechanical mirror and a Fabry-Perot cavity, with k FP given by Eq. 3. k FP scales with L −1 , and so do the other couplings, through their dependence on the cavity volume V c or the single-photon coupling g 0 . We recall that the Fisher information depends on k 2 , which means that it ultimately scales with L −2 . Thus, changing L by an order of 10 will decrease the Fisher information by an order of 10 2 . This contributes to the challenges of realising this scheme. However, it is important to note that there are realistic ways of increasing L without changing the single-photon coupling: One such method was explored in 52 , where L was increased by adding an optical fibre to the cavity.
Furthermore, in the above we proved the optimality of a homodyne detection scheme, but we also found that such a measurement must be performed within a rather narrow temporal window, of timescale 1/k. Let us here estimate how quickly these measurements have to be performed based on the values we calculated for the coupling constant k. The nanospheres displayed the highest single-photon coupling with k Lev ω m ¼ 10 5 Hz for the choice of ω m = 10 2 Hz. Thus any homodyne measurement must be performed within 10 −5 s, so we require at most microsecond precision, which is perfectly achievable. In comparison, we calculated k FP ¼ 2:30 for the levitated microsphere, which allows for a very comfortable ≈0.19 s window.
In spite of these challenges, optomechanical systems come with a number of advantages. They can remain stationary while performing the measurement, in contrast to on-chip BECs or BEC fountains which need to be launched, and the short cycle time of optomechanical systems allows for a large number of measurements to be performed very quickly. An additional point which we did not elaborate on above is that the spatial resolution of optomechanical systems will be extremely high since the oscillator is displaced only by a minuscule distance. As a result, it will be possible to determine very fine local variations in g, something which is not possible using larger systems. The scheme presented in this work also allows for the creation of macroscopic spatial superpositions, which, as pointed out in 11 , is of great interest to testing gravitational collapse models (see for example [53][54][55] ).
Before we conclude, let us now briefly discuss the underlying physical differences between atom interferometry and optomechanical systems for the purpose of gravimetry. We estimate that the QFI for atom interferometry is given by H Q ðTÞ $ n 2 T 4 k 2 c up to an unknown geometric factor, where n is the number of photons that deliver a momentum kick to the atoms, T is the total time over which the gravimetric phase is accumulated, and k c is the laser wavevector. See Supplementary Note 4 for the derivation. If we compare this to the Fisher information for optomechanical systems, we find that the Fabry-Perot cavity has a QFI that is larger by an enhancement factor ξ FP $ c 2 =ðnL 2 ω 2 m Þ, with c being the speed of light. This is due to the cavity confinement, whereby each photon interacts with the oscillator c/ (2Lω m ) times per oscillation cycle, which is also the time period over which the gravimetric phase is accumulated. For the levitated nanosphere, we find a ξ Lev ¼ ξ FP P 2 =ðε 0 V c Þ 2 , where, again, for a micro-object containing~10 13 atoms, the polarizability P is much higher than that of a single atom. In practice, however, both of the enhancement factors will be damped by a factor~1/ (ω m T) 4 with respect to atom interferometry as the time of atomic interferometry T is typically larger than the time 1/ω m of our scheme. Thus the sensitivity Δg is seen to improve by a factor of ffiffiffi n p Lω 3 m T 2 =c $ ffiffiffi n p 10 À4 in our optomechanical scheme with respect to atomic interferometers. As n increases, the differences level out. However, different initial states in an optomechanical system, such as the superposition of two Fock states, will also scale with n 2 (see Supplementary Note 5). Strictly speaking, the enhancement is valid for when the cavity field remains coherent for the time 1/ω m over which our phase accumulation, i.e., κ«ω m (the resolved side-band regime). However, our numerical results indicate that even in the presence of finite decoherence, say, κ~0.1×ω m , the Fisher information is lowered only by a factor of about 10 compared to the case of loss-less cavities. Finally, we can also compare the treatment presented in this work to a position measurement of a classical oscillator that has been displaced due to gravity. While a classical treatment of the problem returns a preliminary measurement sensitivity similar to what we have derived in this work, it fails to take into account effects such as radiation pressure and the full quantum nature of the cavity field. Most importantly, a classical treatment does not utilise the coherent nature of the oscillator, which as we saw above negates any initial thermal noise in the state, and does not allow for the inclusion of other quantum states, such as squeezed states. After completing this work, the authors became aware of similar work carried out by Armata, Latmiral, Plato and Kim 58 .

Methods
Dimensionless Fisher Information. For any numerical calculations, we must ensure that the numerical quantities are dimensionless. In order to calculate I F (t) under decoherence, we separate Eq. 13 above into a dimensionful and dimensionless part by writing Here, ∂ g=∂g ¼ cos θ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m=ð2 hω 3 m Þ p is a dimensionful prefactor. The remaining integral is now dimensionless and is what we will evaluate numerically. A final estimate for Δg can then be obtained by multiplying the value for I F ðtÞ by cos 2 θ m=ð2 hω 3 m Þ, but as this is only a rescaling we chose to only present the results for I F ðtÞ for clarity. Note that for the choice of k ¼ g ¼ 1, this dimensional prefactor is equal to cos 2 θ m=ð2 hω 3 m Þ ¼ 1=g.
Numerical methods. To evolve the system, we use the Python library Qutip 56 and a 4th order Runge-Kutta-Fehlberg method 57 for verification.
If we wish to compute the CFI for states undergoing decoherence, Eq. 15 is no longer valid and we must evolve the state numerically. We do so by computing the dimensionless part I F ðtÞ in Eq. 20 for a mixed state ρ(t). The probability distribution p(x|g) is easy to obtain numerically, since any operator has a finite matrix representation from which we can obtain the eigenstates and use these as our POVM elements. For example, we can define a position operatorx C as a finitedimensional matrix and solve for its eigenstates.
To obtain I F ðtÞ we must also compute the derivative ∂pðxjgÞ=∂ g. This can be done in a number of ways. The simplest one is to use a higher-order method of the central difference theorem. We obtained good and accurate results with the 4th order five-point method. For a function f(x) with parameter x and step-size h, the first derivative with this method is given by f ′ðxÞ ¼ Àf ðx þ 2hÞ þ 8f ðx þ hÞ ½ À8f ðx À hÞ þ f ðx À 2hÞ=ð12hÞ þ Oðh 4 Þ: As this method requires five data point, it is an expensive computation. This was our preferred numerical method as computing the CFI can still be done within reasonable timescales using the optimised master equation solver provided by the Qutip library. It does however contain two different sources of numerical errors: errors in the solver and errors that originate from the cut-off in the numerical derivative.
To verify that the error introduced by the numerical differentiation is not severely affecting the results, we used an additional method which provides an exact result. We reproduce it here for completion and in the hope that it might benefit others. We start by noting that as long as the POVM element П x does not depend on g, we can write the derivative as Note that we are differentiating with respect to g instead of g and that we have suppressed the dependence of t for clarity. This statement also holds for subsystems of ρ(g), which we can see by noting that the derivative distributes over a joint separable system ρ AB ¼ ρ A ρ B as Performing a measurement with П x that only acts on subsystem A then gives The second term reduces to zero because tr½∂ g ρ B ¼ ∂ g tr ½ρ B ¼ 0. While we have shown this for separable states, the same argument can be extended to entangled states by linearity.
In order to obtain the evolution for this state, we must now solve a modified version of the master equation. That is, given the Lindblad equation in Eq. 19, _ ρð gÞ ¼ À i hĤ ð gÞ; ρð gÞ Â Ã þLρð gÞL y À 1 2 fρð gÞ;L yL g; ð27Þ where fÁ; Ág denotes the anti-commutator, we now differentiate both sides with respect to g to obtain ∂ g _ ρð gÞ ¼ À i h ∂ gĤ ð gÞ; ρð gÞ h i À i hĤ ð gÞ; ∂ g ρð gÞ h i À 1 2 f∂ g ρð gÞ;L yL g; where we have again used the notation ∂ g ¼ ∂=∂ g. A more complicated form is obtained if the Lindblad operatorsL depend on g, which here is not the case. In coupled form, we can write þL∂ g ρL y À 1 2 f∂ g ; ρL yL g 0 @ 1 A : The system can be solved using any standard higher-order method, such as the family of Runge-Kutta ODE solvers. Note that the Qutip Master Equation solver cannot be used as Eq. 27 is not in standard Hamiltonian form.
Once the time-evolved state ∂ρ=∂ g has been obtained, we proceed as usual to compute the probability distribution with the set of POVM elements {П x }. With this method, we avoid round-off errors that appear in the five-point numerical derivative above.
Numerical stability. Let us make a few remarks regarding the numerical stability of the simulation. We start by considering the nature of coherent states and how they are represented numerically. Coherent states have support on infinite Hilbert spaces, whereas numerically we must work with finite matrices. It is therefore necessary to introduce a cut-off in the dimension used to represent the state. This leads to a gradual loss of coherence as information is pushed beyond the cut-off. In other words, numerically we use a finite Hilbert space H, meaning that we truncate the space by lettingâ y jN À 1i ¼ 0, where dimðHÞ ¼ N. Furthermore, the appearance of ðâ yâ Þ 2 inÛðtÞ causes the system to become anharmonic and numerical instabilities grow fast for Hilbert spaces with small dimension N < 50.
The amount of information lost when using smaller Hilbert spaces is difficult to assess, since any good ODE solver will preserve the purity of the state throughout the simulation. Rather, the loss of information can be noted as a gradual deterioration of the trajectory in phase space, with the effect that states fail to return to their original position in phase space at t = 2π. That is, we require that hxð0Þi % hxð2πÞi and hpð0Þi % hpð2πÞi for the simulation to be deemed stable. This is however only true for noiseless evolution.
The system dynamics depend strongly on the dimensionless constants k and g. Larger k and g will cause the system to evolve more rapidly, as evident from their appearance in the phase of the state in Eq. 7. This in turn causes the numerical inaccuracies to accumulate more rapidly. When computing the CFI for mixed states, we restrict our investigations to the parameter range k ¼ g ¼ 1 for precisely this reason.
Finally, it should be noted that we have not provided a full error estimate for any of the results computed here. However, since we are only interested in the general behaviour of the CFI, any small inaccuracies in the numerical estimates will not matter for the results obtained in this work.
Data availability. The datasets generated during and/or analysed during the current study are available in the Coherent-states-Fisher-information repository on Github at (https://github.com/sqvarfort/Coherent-states-Fisher-information).