Constraining modified gravity with quantum optomechanics

We derive the best possible bounds that can be placed on Yukawa- and chameleon-like modifications to the Newtonian gravitational potential with a cavity optomechanical quantum sensor. By modelling the effects on an oscillating source-sphere on the optomechanical system from first-principles, we derive the fundamental sensitivity with which these modifications can be detected in the absence of environmental noise. In particular, we take into account the large size of the optomechanical probe compared with the range of the fifth forces that we wish to probe and quantify the resulting screening effect when both the source and probe are spherical. Our results show that optomechanical systems in high vacuum could, in principle, further constrain the parameters of chameleon-like modifications to Newtonian gravity.


Introduction
General Relativity is one of the most successful theories of nature, but there are compelling reasons to explore modifications to the behaviour of gravity on both large and small scales. Most of the precise predictions of General Relativity have consistently been demonstrated experimentally: among many others these include the perihelion shift of Mercury [1] and the existence of gravitational waves [2]. Similarly, the current standard cosmological model, the Λ Cold Dark Matter (ΛCDM) model, is another of General Relativity's success stories. However, in order to match observation, ΛCDM requires a positive cosmological constant [3,4]. This is backed up by observations of supernovae, which indicate that the Universe's expansion is accelerating [5]. While a natural part of General Relativity, a cosmological constant poses a theoretical challenge to particle physics since the small observed value is inherently sensitive to high-energies, requiring delicate balancing [6]. Furthermore, many theories of high energy physics that attempt to solve this and other problems -such as building a consistent quantum theory of gravity -predict deviations from General Relativity. These theories are collectively known as modified gravity theories.
Modified gravity theories, however, typically face a difficult challenge in the form of solar system tests of Newton's laws. Models that differ from General Relativity significantly enough to explain the observed acceleration of the Universe on large scales are typically ruled out by their predicted deviations on smaller scales (solar system and laboratory tests) [7][8][9]. There are a large variety of approaches to modified gravity -see Koyama [10] for a comprehensive review -but many models attempt to address the problem of solar system tests via a screening mechanism [11]. Such mechanisms can be built into modified gravity theories to conceal deviations on solar system scales, without changing the large scale behaviour. An approach considered by many authors is the chameleon mechanism [12][13][14]; the basic idea is to add a scalar field that couples directly to gravity in a manner that depends on the local density of matter. In high-density regions, such as inside a galaxy, the effects of modified gravity are screened out, allowing the theory to evade solar system tests. In the low-density void regions between galaxies, however, the effects of modified gravity would be unscreened.
If such a density-dependent gravity mechanism is at play, it ought to be detectable in principle by high-precision laboratory experiments. In particular, the fundamental sensitivity improvements offered by quantum systems are especially promising [15]. At the moment, the detection of modified gravity, and in particular, chameleon fields, has been explored through a diverse variety of methods. Searches with classical systems include theoretical proposals for torsion balance tests of fifth forces [16][17][18][19][20][21][22], some of which have already been carried out as experiments [23,24]. Additional proposals suggest that experiments which measure Casimir forces may also be used to constrain chameleon theories [18,[25][26][27][28]. In atom interferometry, which is already routinely used for quantum sensing, the uniformity of the atoms as well as the additional sensitivity gained from the superposition of flight-paths has led to impressive precision gravimetry sensitivities [29][30][31][32]. Several proposals have explored in depth the possibilities of searching for modified gravity and dark energy with atom interferometry [33][34][35][36][37][38][39], and some of the most stringent bounds on existing theories have been obtained in this way [40,41]. Further viable routes towards detecting modified gravity include ultra-cold neutron experiments [42][43][44][45][46][47][48] and neutron interferometry [46,[49][50][51][52]. Finally, tests of atomic transition frequencies [53,54], close examination of vacuum chambers and photo-detectors [55,56], as well as tests of the electron magnetic moment [57] have also been proposed.
An additional approach to detecting the small-scale effects of modified gravity and screening is to take advantage of recent developments in the field of optomechanics, where a small mechanical element is coupled to a laser through radiation-pressure [58,59]. Optomechanical system encompass a diverse set of platforms which range from microscopic movable mirrors as part of a Fabry-Pérot cavity [60], levitated particles [61], clamped membranes [62], liquid Helium [63] and trapped cold atoms [64]. When the mechanical element is cooled down to sufficiently low temperatures, it enters into a quantum state that can be manipulated through measurements and optical control techniques. Ground-state cooling has been demonstrated across a number of platforms, including clamped membranes [65,66] and recently also for levitated systems [67]. Optomechanical systems show promising potential as both classical and quantum-limited sensors [68][69][70], and recent studies have proposed their use as gravity sensors [71][72][73][74]. In fact, experimental searches for fifth forces with classical optomechanical setups have already been performed (see e.g. [75,76]), where the bounds achieved fell within those excluded by atom interferometry. A key question, which we explore in this work, therefore becomes whether an optomechanical sensor in the quantum regime can improve on these Vacuum chamber Piezo Figure 1. A gold source mass attached to a shear piezo oscillates to create a time-varying gravitational field. The field, which potentially contains deviations from Newtonian gravity, is detected by an optomechanical probe system where the photon numberâ †â couples to the mechanical positionx mech asâ †âx mech , here presented as a moving-end mirror in a Fabry-Pérot cavity. The amplitude x 0 of the source mass oscillation is a fraction of the total distance x 0 between the systems. By accounting for the vacuum background density, we may also compute bounds on the parameters of the chameleon screening mechanism.
bounds. For an overview of searches for new physics with levitated optomechanical systems, see the recent review by Moore et al. [77]. The advantage of optomechanical sensors, as opposed to, for example, cold atom interferometry is that the sensitivity of the system can be improved while retaining the compact setup of the experiment. In contrast, improving the sensitivity of atom interferometry primarily relies on increasing the length of the flight-path of the atoms.
The key question we seek to answer in this work is: what fundamental range of parameters of modified gravity theories could ideally be excluded with a quantum optomechanical sensor? To address this question, we consider an idealised system described by a nonlinear, dispersive, optomechanical Hamiltonian which couples the optical and mechanical degrees of freedom through a nonlinear radiation-pressure term. This Hamiltonian is often linearised for a strong coherent input drive, however the fully nonlinear (in the sense of the equations of motion) Hamiltonian is a more fundamental description. While all quantum systems are affected by noise, we here assume that the coherence times can be made long enough for the measurement protocol to be carried out. As a result, our analysis explores the bounds in the absence of environmental noise and decoherence. We then consider the gravitational field that arises when a source mass is placed next to the sensor.
Since it is often difficult in experiments to distinguish a signal against a constant noise floor, we consider an oscillating source mass, which gives rise to a time-dependent gravitational field. Such a signal can then be isolated from other common low-frequency 1/f noise sources via a Fourier analysis of the data. To determine whether our analysis is valid in the case of a chameleon field, we derive the time-dependent potential that results from the source mass from first principles, where we find that a potential that moves with the mass is the correct choice for non-relativistic velocities. Another key consideration for optomechanical systems is the relatively large size of the optomechanical probe. This has been found to be significant in previous classical experiments with chameleon fields, such as the MICROSCOPE experiment [78,79], and we find that it also contributes significantly to the chameleon screening of the fifth force in the envisioned setup of the quantum experiment we consider here (as opposed to, for example, cold atoms, where the screening length of the atomic probes is very small). To take the finite screening length into account, we go beyond the common approximation that the probe radius is small compared to the range of the chameleon field and derive analytic expressions for the modified force seen by the probe.
Then, using tools from quantum information theory and quantum metrology such as the quantum Fisher information, we are able to estimate the fundamental sensitivity for detecting deviations from Newtonian gravity. To further improve the sensitivity, we also consider known ways to enhance the optomechanical sensor in the form of squeezed light and a modulated optomechanical coupling [74]. Our main results include the bounds presented in figure 4, which shows the parameter ranges of modified gravity theories that could potentially be excluded with an ideal optomechanical sensor. The bounds are computed for a specific set of experimental parameters. To facilitate investigations into additional parameter regimes, we have made the code used to compute the bounds available (see the Data Availability Statement). While experiments are unlikely to achieve the predicted sensitivities due to noise and systematic effects, our bounds constitute a fundamental limit for excluding effects beyond Newtonian gravity given the experimental parameters in question.
This work is structured as follows. In section 2 we present the proposed experimental setup and optomechanical Hamiltonian, and then we proceed to discuss Yukawa potentials as a modification to the Newtonian gravitational potential in section 3. We consider those sourced by a chameleon field and provide a first-principles' derivation of the time-dependent potential that results from the mass oscillating around an equilibrium position. We also discuss screening effects inherent to chameleon fields and derive the screening effect that arises from the size of the optomechanical probe. In section 4, we linearise the modified gravitational potential, and in section 5, we provide an introduction to quantum metrology and the quantum Fisher information. These tools allow us to present analytic expressions for the fundamental sensitivity of the system, which we do in section 6. The work is concluded by a discussion in section 7 and some final remarks in section 8.

Optomechanical model and dynamics
In this section, we introduce the model of the optomechanical system and show how the effects of a timevarying gravitational field can be included in the dynamics.

Experimental setup
We envision an experimental setup similar to that used in [80], where an oscillating source mass made of solid gold is placed in a vacuum chamber adjacent to an optomechanical probe (see figure 1). We have chosen gold because we require the highest possible density in order to detect gravitational effects and maximise the effect of density-dependent screening mechanisms such as chameleon fields ‡. The source mass oscillates back and forth, which can be achieved in a number of different ways [81]. One such implementation is with the help of a shear piezo, which oscillates at a fixed frequency. The optomechanical probe is then allowed to move along the same axis as the oscillating mass. By injecting light into the cavity, the position of the optomechanical coupling is dispersively coupled to the optical field through radiation-pressure. The light then picks up a phase shift conditioned on the displacement of the mechanical mode, which has been influenced by the gravitational force. Therefore, information about the gravitational field is imprinted on the optical state. The light is then collected and measured either as it leaks from the cavity or through a scheme where the cavity is coherently opened to access the full intra-cavity state [82].
While the optomechanical interaction can generally be described with the same dynamics for a large range of systems, the force and strength of the coupling differ for each platform. In this work, we begin with a general description of a single interacting mode, but later specialise towards a spherical mechanical element since it allows for analytical treatments of some modified gravity potentials.
The optomechanical Hamiltonian, which governs the dynamics of the optomechanical probe, is given by (in the absence of an external gravitational field): where ω c and ω mech are the oscillation frequencies of the optical cavity mode and mechanical mode respectively, with annihilation and creation operatorsâ,â † andb,b † . We have also definedN a =â †â and N b =b †b as the photon and phonon number operators.
The coupling k(t) is the (potentially time-dependent) characteristic single-photon interaction strength between the number of photons and the position of the mechanical element. It takes on different forms depending on the optomechanical platform in question. Among the simplest couplings is that for a moving mirror, of mass m, that makes up one end of a cavity, k = (ω l /L) /(2mω mech ), where ω l is the laser frequency and L is the length of the cavity. In this work, we also consider modulating the coupling in time; it has been previously found that such modulations can be used to enhance the sensitivity of the optomechanical sensor if they can be made to match the oscillation of the external force [74]. Modulation of the optomechanical coupling can be introduced in different ways depending on the experimental platform in question. For example, the mechanical frequency of a cantilever can be modified by applying an oscillating electric field [83,84], and a modulated coupling arises naturally through the micro-motion of a levitated system in a hybrid electro-optical trap [85][86][87].
All quantum systems are affected by noise due to their interaction with the environment. Such an interaction usually results in dissipation and thermalisation, which in turn leads to decoherence of the offdiagonal elements of the quantum state. For cavity optomechanical systems, common sources of noise include photons leaking from the cavity, as well as thermlisation of the mechanical element due to interactions with the surrounding residual gas, or from vibrations from the mount [59]. The nature of the noise is unique to each experimental platform and must be carefully modelled in each case.
In this work, we are interested in deriving the best-possible sensitivity that an optomechanical system can achieve. We therefore assume that the Q-factor of the cavity is high enough that the system stays coherent throughout the duration of our measurement protocols. Recently, Q-factors of 10 9 have been demonstrated in magnetically levitated meso-mechanical systems [88], and linewidths of 81 ± 23 µHz have been measured [89]. We also assume that the system has been cooled to temperatures such that the surrounding environment does not cause the mechanical mode to heat up during the protocol. To reduce unwanted vibrations or gravitational noise, it is also possible to add decoupling stages in the experiments [90], such as suspension stages made by fused silica fibres [91,92]. Under these conditions, it is possible to consider an approximately unitary description of the experiment, which we shall use to derive a fundamental limit of the sensitivity that could in principle be achieved with an optomechanical system. To then describe a realistic experiment, all of the above effects must be taken into account. We discuss this and other potential future work in section 7. When treating the system in a closed and ideal setting, we can model the initial state as a separable state of the light and the mechanical element. For the optical state, we consider injecting squeezed light into the cavity. Squeezed light has been shown to fundamentally enhance the sensitivity to displacements [15]. By including squeezing here, we generalise our scheme to include these input states. However, we note that in order to improve the sensitivity overall, it is always more beneficial to increase the number of photons rather than squeezing the system. Squeezing also reduces quadrature noise [93]. The state of the mechanical element, on the other hand, is most accurately described as thermal at a non-zero temperature. With these assumptions, the initial state of the system can be written aŝ where |ζ =Ŝ ζ |µ c is a squeezed coherent state of the optical field whereŜ ζ = exp (ζ * â2 − ζâ †2 )/2 and where the coherent state satisfiesâ |µ c = µ c |µ c . The squeezing parameter can also be in spherical polar form as ζ = r sq e iϕ . Squeezed states can be generated through four-wave mixing in an optical cavity [94] or parametric down-conversion [95]. See also Ref [96] for a review of squeezed state generation. The parameter r T of the thermal state arises from the Bose-Einstein distribution and is defined by tanh r T = exp − ω mech 2 kB T , where T is the temperature of the system and k B is Boltzmann's constant.

Modelling the gravitational force
In order to compute the sensitivity bounds for detecting modified gravity, we model the effect of the gravitational force from the moving source mass on the optomechanical system as a contribution to the dynamics. When the force is weak, it can be linearised and included as a displacement term in the optomechanical Hamiltonian in equation (1). In this section, we provide a general derivation of this linearised force, while in section 4 we specialise to Yukawa-like and chameleon modifications to the Newtonian potential.
The linearisation is necessary to properly describe the quantum dynamics of the setup with the current theoretical machinery; we will however describe the chameleon field in full generality to allow for future work to improve the theoretical description.
We start by assuming that the source mass and the mechanical element of the optomechanical system are constrained to move along the x-axis. We let the mechanical element be subject to a harmonic potential centered at x = 0 and we label the position of the source mass x S (t). Then, we assume that there is a small perturbation to the centre-of-mass position of the mechanical element that we call δx, and, assuming that |x S (t)| |δx| at all times, we write the relative distance between the systems as x S (t) − δx. Provided that δx remains small, we can Taylor expand the system to first order in δx. Given a generic potential term V (x S (t) − δx), we find, to first order in δx: The first term in equation (3) represents a time-dependent shift of the overall energy, however it does not depend on the position of the optomechanical system. The second term describes the (potentially timedependent) displacement of the mechanical element with δx. The second-order term in δx leads to a shift in the mechanical frequency that we do not model here, but dynamics of this kind have been previously studied [97]. The expansion in equation (3) is valid as long as δx remains small such that the higher-order terms can be neglected. We outline the conditions for this being true in the Discussion (see section 7). We proceed to promote δx to an operator δx →x mech , which can be written in terms of the annihilation and creation operatorsb andb † of the mechanical element aŝ where x zpf = /(2 mω mech ) is the zero-point fluctuation of the mechanical oscillator. It should be noted here that the dynamics of a nonlinear optomechanical system with a driving term proportional tox 2 mech has been solved, however the inclusion of these effects adds significant complexity the mathematical treatment of the system [97], while it will likely not result in a significant improvement of the sensitivity.
The full optomechanical Hamiltonian including the modified gravitational potential can then be written asĤ whereĤ 0 is given in equation (1) and where the time-dependent modified Newtonian gravitational force is contained in the second term. The time-evolution of the system with the Hamiltonian in equation (5) can be written as the following time-ordered exponential:Û where the time-dependence of the gravitational potential inĤ(t) requires careful consideration. Such dynamics have been studied previously [97,98] and provide a short overview of the treatment Appendix C. Later on in this work, we use the expression forÛ (t) to derive the sensitivity of the system to modifications of the Newtonian potential, but first, we will study the form of the modifications in-depth.

Modified gravitational potential and screening from the source and the probe
In this section, we discuss an example of how the chameleon mechanism would alter the Newtonian force law on sub-millimetre scales. We write all equations in terms of SI units, but energy units (as used elsewhere in the literature) can be restored by setting = c = 1 throughout.

Yukawa modifications to the gravitational force law
Although there are many ways of modifying Newton's laws on short distances, perhaps one of the best motivated theoretically is to add a Yukawa term to the potential. Yukawa potentials are ubiquitous in scalar field theories, since they are the solution to the sourced (inhomogeneous) Klein-Gordon equation for a massive field in the case of spherical symmetry. As a consequence of Lovelock's theorem [99], modifications to General Relativity require either additional degrees of freedom such as a scalar field, or more exotic scenarios such as large extra dimensions, higher derivatives, or non-locality. Consequently, additional scalar fields are common in modified gravity theories. These act like a fifth force, and for a source of mass M S and test particle mass m, give rise to a gravitational potential of the form: where α parametrises the intrinsic difference in strength between the Yukawa-like fifth force and gravity, while λ parametrises the range of this fifth-force. Note that it is possible to have α 1 and still agree with existing constraints, provided the force is sufficiently short range to have evaded tests of gravity on short distances.
For this work, we will consider a chameleon screening mechanism that gives rise to a Yukawa-like force. However, the methods we describe can be broadly applied to many different Yukawa-type modifications of the gravitational field on short distances. In the chameleon mechanism, short distance modifications to the Newtonian force law are screened from the reach of solar system tests by the presence of a density-dependent scalar field, known as the chameleon field. In regions of relatively high average density -such as can be found inside a galaxy -the chameleon field has a high mass, making it hard to detect at colliders and altering the gravitational force law in such a way as to be consistent with solar-system experiments (this is the 'screening' effect). However, in regions of low-density -such as in cosmological voids -the field is lighter and the effects of modified gravity unscreened. This allows modified gravity theories to have substantial effects on cosmological scales, while being difficult to detect on galactic or solar-system scales.
We review the properties of chameleon fields in Appendix A. The net effect of the chameleon scalar field φ is to modify the effective Newtonian potential affecting a test particle. Specifically, the effective potential at position X is given by where Φ N is the standard Newtonian potential, and Φ C is the modification to it arising from the chameleon field. The parameter M (here chosen to be a mass to give the correct units for a potential) determines how strongly the chameleon field affects test particles and arises from the non-minimal coupling of the chameleon field to curvature as discussed in Appendix A.
In this work, we consider a chameleon model with an effective interaction potential We explore only the case n = 1 in this work: other models and choices of n are possible, but we choose this specific example to demonstrate how the method works in principle. This model has two parameters; Λ, which characterises the energy scale of the chameleon's self-interaction potential; and M which is defined above. For n = 1 the background value of the field, φ bg , in an environment of constant mass density ρ bg is given by In the centre of the source, the chameleon field reaches its minimum value of φ S (which can be obtained by replacing the density ρ bg in equation (10) with the source density ρ S ). The mass of the chameleon field, m bg , is density dependent (see Appendix A) and given by The key question for us is how the field results in a force on the optomechanical sensor. This is what we consider next.

Force on the optomechanical sensor
The effect of a chameleon field is in principle detectable in a high-vacuum environment. In practise, this requires extremely precise acceleration measurements, which an optomechanical system can provide. While the optomechanical probe can come in many diffferent shapes, in this work, for the sake of simplicity, we model both the source mass and the detector probe as spheres. This allows us to compute the sensitivity using the chameleon force between these two spheres. There are therefore two effects to consider: the response of the field in equation (8) to the spherical source, and the response of the probe to that field. Due to the nature of the chameleon field, a non-point-like probe will not simply follow the gradient of equation (8) as would a test particle: instead there is an additional screening effect due to the interactions of the probe itself with the field.
To derive the force that acts on the sensor, we consider the field inside the vacuum chamber. Burrage et al. [33] derived the chameleon field around a spherical source of mass M S and radius R S as a function of distance from the centre of the sphere, r. They assumed in their derivation that the range of the chameleon force was large compared to the size of the source (that is, m bg R S c/ 1). To allow us to consider a broad parameter space, we do not assume that either m bg R S c/ 1 or m bg R P c/ 1 where the indices S and P denote the source and probe, respectively. In what follows, we go beyond existing studies in this regard by including sources (probes) with non-negligible size compared to the range of the force.
We use the same asymptotic matching approach as Burrage et al. [33] to obtain an expression for the chameleon field around a spherical matter distribution: Here, i = S, P for the source and probe respectively. φ i is the equilibrium value of the chameleon field in a material with the source (probe) density: the field will attain this value at some radius S i ≤ R i . There is then a transition layer where the field φ i increases to its surface value, before increasing to the equilibrium value in the background density, φ bg . The range of the force outside the source is controlled by the density-dependent chameleon mass, m bg . The full derivation of equation (12) is given in Appendix A.1.
The value of the length scale S i depends on the source/probe properties, the chameleon model, and environmental properties. For the model we consider here, it is found by solving the following cubic equation: In the m bg R i c/ → 0 and φ bg φ i limits, this reduces to which is the result found by Burrage et al. [33]. S i parametrises the screening effect of the chameleon mechanism for a spherical source/probe: for example, when S i is much lower than R i , the field is effectively unscreened while for S i ≈ R i the field is heavily screened. Outside of the source (r > R i ), we see from equation (12) that the scalar field, and thus the modified gravitational potential, has an effective Yukawa form. Thus, a chameleon field of this type would manifest as a Yukawa-like modification to the acceleration of a test particle. Since our proposed experimental setup will involve measuring acceleration outside of the sphere, we only need the r > R i part of the solution.
When viewed as a Yukawa-type force of the form considered in equation (7), and when the probe itself does not contribute to the screening, the resulting fifth-force strength α and range λ are given by where M P ≈ c/(8πG) = 4.341 × 10 −9 kg is the reduced Planck mass (here expressed as a mass rather than an energy), α bg depends on the background density through ξ S , which is given by [33] As long as the optomechanical sensor is approximated as a point-particle, such that R P /λ bg 1, the force felt by the optomechanical probe can therefore be written as where X S is the vector-position of the source. The point-particle approximation is however quite severe, especially for an optomechanical probe, the radius of which can be quite large compared with the range of the force in question. We proceed to consider the screening from the probe in the following section.

Chameleon screening from the optomechanical probe
Compared with the atoms used in alternative approaches to the detection of fifth-force modifications to gravity, such as atom interferometry, the optomechanical probe can potentially be relatively large compared to the range of fifth forces. This can result in significant contributions to the chameleon field screening. The screening depends strongly on the geometry of the system; in general, numerical methods are needed to compute the full screening [100]. As such, it is difficult to estimate the screening for say a Fabry-Pérot moving-end mirror; however the problem is simplified when both the source sphere and the probe are spherically symmetric. This is the case when the mechanical element in the optomechanical system is a levitated sphere, made, for example, by silica.
To estimate the extent of the screening for a spherical optomechanical probe, we consider the force that arises from the movement on the time-dependent mass. See Appendix B.3 for the full calculation. In the limit where the probe radius is much smaller than the distance between the probe and the source sphere R P |X S (t)|, we find the following expression for the force: where the sensor-dependent fifth-force strength is defined as where we have added the subscript 'P ' to denote that screening from the probe is here taken into account. Furthermore, ξ S and ξ P (again labelled S for the source and P probe, respectively), are given in equation (16).
To compute ξ P , we replace M S , R S and ρ S with M P , R P and ρ P . Finally, the function f is a form-factor given by This approaches 1 in the x = m bg R P c/ = R P /λ bg → 0 limit, in which case equation (18) reduces to the result of Burrage et al. [33] for the force between two spheres. Since spherical probes or source masses generally maximise the screening [100], equation (18) can be interpreted as a conservative estimate of the screening due to the shielding from the probe. Burrage et al. [33] make use of the R P /λ bg → 0, since the probe radius in the case of atom interferometry is typically much smaller than λ bg . For an optomechanical probe, however, the additional screening introduced by the probe can be substantial, but not, as we shall see, detrimental. In what follows, we compute the sensitivity both with and without the screening from the probe, where the latter corresponds to approximating the probe as a point particle.

Potential from a moving source-mass
In this work, we consider a moving source mass. This brings up a consideration of how the chameleon field responds to the motion of the source mass. For the gravitational field, we know that changes in the potential propagate outwards at the speed of light, and thus the appropriate potential to use is the retarded Newtonian potential. The situation is less clear for the scalar field, however. Since it is massive, it is not immediately obvious information will propagate outwards at the speed of light. To get an idea of its behaviour we need to know the speed, v I , at which information propagates through the scalar field. We show in Appendix B.2 that this is also, in fact, the speed of light. Consequently, the potential at 3D position X can be approximated by the time dependent form where t should be replaced with the retarded time given by equation (B.9), however, we can ignore this for the non-relativistic speeds and distances considered in this setup. Since both the chameleon field, φ, and the metric, g µν (which gives rise to the Newtonian potential, Φ N ) are well-defined dynamical quantities, the time-dependence of this potential is well-defined. We note at this point that if quantum corrections are large, the effective speed of information propagation for the scalar field, v I , may differ from c [101]. However, large quantum corrections of this size would mean that we cannot readily use the effective field theory treatment of the chameleon field assumed throughout [102], so we do not consider this effect here. We can therefore use equation (21) in the discussion that follows to measure the values of α and λ, and thus the parameters Λ and M of the chameleon field.

Linearised modified Newtonian potential
In order to compute the sensitivity of the optomechanical system, we need to include the force on the sensor shown in equation (18) into the dynamics of the optomechanical system. It is possible to obtain the solution numerically, but in order to obtain analytic expressions, we choose to linearise the Yukawa modification of the force for small oscillations of the source-mass. We let the time-dependent distance between the systems x S (t) be given by: where is a dimensionless oscillation amplitude defined as a fraction of x 0 , where ω 0 is the oscillation frequency and φ 0 is a phase shift that we specify later in order to maximize the sensitivity. In the following two sections, we show the linearisation of the force for a generic Yukawa potential, and for the chameleon force with a large optomechanical probe that contributes to the screening.

Linearising the Yukawa potential
We now linearise the contributions from the Yukawa potential to equation (7) for small oscillation amplitudes 1. We note that, for specific values of α and λ, higher order contributions to the Newtonian gravitational force may be larger than the first-order contributions to the Yukawa force. It is therefore important that, when taking data in an experiment, we determine the origin of the observed values (see the Discussion in section 7). Linearising, we obtain where g N = G M S /x 2 0 is the Newtonian gravitational acceleration at the equilibrium distance x 0 and where we defined the two parameters which quantify the deviation of the constant and the time-dependent part of the force from the Newtonian one, respectively. We make this distinction because in an experiment, it is often possible to isolate a timedependent signal from a constant noise floor. In addition, systematic effects such as the Casimir effect can be effectively screened out in this way (we return to this point in the Discussion in Section 7). We will therefore focus on estimating σ as part of our analysis.

Linearising the chameleon potential from a screened spherical probe
For a spherical optomechanical probe, the force between the probe and the source is given in equation (18). We note that the form factor shown in equation (20), which arises due to the screening from the optomechanical probe, depends on x S (t) and is therefore time-dependent. To determine the constant and time-dependent contributions, we again assume that the source sphere oscillates around the equilibrium distance x 0 according to equation (22). Noting that x S (t) only enters into the y-terms in equation (20), we find where the expressions for κ and σ now read where we have defined: The expressions in equation (27) arise from the form-factor in equation (20). For the parameter regimes considered in this work, we find that R P /λ bg 1. This means that the form factors become A(R P /λ bg ) = 1 and B(R P /λ bg ) = 0. As a result, κ and σ simplify to which has the same form as equation (24). We are now ready to compute the sensitivities of the optomechanical sensor, but first, we provide a brief introduction to the quantum metrology tools we use for this purpose.

Quantum metrology and ideal bounds
In this work, we are interested in the best-possible sensitivity that can be achieved with the optomechanical probe. To determine the sensitivity of the probe, we turn to tools from quantum metrology. Specifically, we focus on computing the quantum Fisher information (QFI), which we denote I θ , where θ is the parameter that we wish to estimate. Intuitively the QFI can be seen as a measure of how much the quantum state of the system changes given a specific encoding of θ. The QFI then provides a measure of the change in the state with θ compared with the case when the state is unaffected. See also Ref [103] for an intuitive introduction to the QFI and related concepts in quantum metrology. The connection to sensitivity stems from the fact that the QFI provides a lower bound to the variance Var(θ) of θ through the quantum Cramér-Rao bound [104,105]: where M is the number of measurements or probes used in parallel. The standard deviation of θ is then given by ∆θ = 1/ √ M I θ .
For unitary dynamics and mixed initial states written in the form ofˆ = n λ n |λ n λ n |, the QFI can be cast as [106,107]: where the operatorĤ θ is defined asĤ θ = −iÛ † θ ∂ θÛθ . Here,Û θ is the unitary operator that encodes the parameter θ into the system.
In our case,Û (θ) is the unitary operator that arises from the Hamiltonian in equation (5), and the effect we wish to estimate is the effect of the Yukawa potential on the probe. Therefore, in order to compute I θ , we must first solve the time-evolution of the system, which is often challenging when the signal is time-dependent, as is the case for us here. Some of these challenges can however be addressed by making use of a previously established method for solving the Schrödinger equation using a Lie algebra approach [108]. Details of this solution were first used to study a purely Newtonian time-dependent gravitational potential [74] and can be found in Appendix C.
Using the expression for I θ in equation (30), we can derive a compact expression for the QFI that represent the sensitivity with which modifications to Newtonian gravity can be detected. In our case, we let the parameter θ of interest be either κ or σ as defined in equation (24). By then applying the Cramér-Rao bound, we can derive the standard deviation for each parameter. We then consider the ratios ∆κ/κ or ∆σ/σ, which describe the relative error of the collective measurements.
In this work, we say that we can distinguish modifications to the Newtonian potential if the error in κ and σ is smaller than one, that is, when ∆κ/κ < 1 or ∆σ/σ < 1. Note that, to find the sensitivity to the actual values of, for example, α and λ, we would need a full multi-parameter likelihood analysis, which requires us to go beyond the regular error-propagation formula for the parameter we consider here. Such an analysis is currently beyond the scope of this work. Instead, we focus mainly on detecting σ, since it is the amplitude of the time-dependent signal.
Unfortunately, the QFI does not actually reveal the optimal measurement that saturates the quantum Cramér-Rao bound. To obtain this information, one must compute the classical Fisher information for a particular measurement and examine whether it saturates the quantum Fisher information. It is known that, when the optomechanical coupling is constant and takes on specific values, that a homodyne measurement of the optical field is optimal [71,74]. When the optomechanical coupling is modulated at resonance, as is the case here, the optimal measurement is not yet known. The gravitational interaction between the source and the optomechanical probe results in a phase shift of the optical state. Therefore, the utility of a homodyne measurements can be expected also for the case of modulated optomechanical coupling, but we leave this specific analysis to future work. In practise, once the optomechanical probe has interacted with the source, the system is measured to extract information about the gravitational force. Standard measurements that are performed on the optomechanical system include homodyne and heterodyne measurements of the cavity field, as well as photon detection measurements, which can either be resolving (counting the number of photons) or non-resolving (merely detecting the presence of a photon). In a homodyne measurement, the output light from the optomechanical system is brought into interference with a local oscillator light field which comes from the same source as the input light field of the optomechanical system. This is the same measurement principle that is, for example, employed in a Mach-Zehnder interferometer to infer a phase shift on a light field. Heterodyne measurements, on the other hand, compare the collected light with a different coherent state reference. The usefulness of each measurement depends on the situation at hand. Since we focus on deriving the best-possible sensitivities in this work, we leave it to future work to analyse the sensitivity that can be gained from specific measurements.

Results
We are now ready to compute the sensitivities that can be achieved with an ideal optomechanical sensor for detecting modifications of gravity. Specifically, we consider a region of parameter space to be possible to exclude using the optomechanical sensor when the best precision possible on the parameters α and λ (or the chameleon parameters M and Λ) is sufficient to distinguish them from zero, their values in ordinary General Relativity. The filled-in contours show the radio without screening from the probe. The lines instead show the ratio for when the probe is spherical and contributes to the screening. As a result, the strength of the force is reduced. The parameters used to make these plots can be found in table 1. The mapping between the (α, λ) and (M, Λ) spaces is non-trivial. The left hand figure shows the range of force modifications that can potentially detected, which is different to the range of theoretically-interesting chameleon parameters, shown on the right.

Fundamental sensitivities
We first present some simple expressions for the sensitivities that can be achieved, and we then proceed to compute the parameter regions that could potentially be excluded with an optomechanical sensor. When the source mass oscillates at the same frequency as the optomechanical system, that is, when ω 0 = ω mwch , the effects accumulate and cause the position of the optomechanical system to become increasingly displaced.
Following the outline in Appendix C we find the following expressions for the sensitivities for κ and σ at time t ω mech = 2πn (see [74] for a detailed derivation). For large enough temperatures in the mechanical state, such that r T 1, the expressions simplify and we find that the sensitivities ∆κ and ∆σ are given by where n is an integer, and for an optomechanical coupling k(t) ≡ k 0 and phase φ 0 = π, and where the variance (∆N a ) 2 of the photon number is given by [74] ( where r sq and ϕ are the squeezing amplitude and phase, and where µ c is the coherent state amplitude of the optical mode. The expression in equation (33) is maximised when e −iϕ/2 µ c is completely imaginary, which causes the last term of equation (33) to vanish. This can be achieved by assuming that µ c ∈ R and setting the squeezing phase to ϕ = π/2. The other parameters in equation (33) have been previously defined in the text (see also table 1 for a summary). The sensitivities can be improved by modulating the optomechanical coupling at the same frequency as the gravitational signal [74]. In this work, we choose a sinusoidal modulation with k(t) = k 0 cos(ω k t), where k 0 is the amplitude of the modulation and ω k is the modulation frequency. At resonance, when ω k = ω mech , and for the optimal phase choice φ 0 = π/2, we find that the sensitivities for measuring κ and σ become Here, equation (35) scales with n −2 rather than n −1 . This enhancement arises from the additional modulation of the optomechanical coupling, and was already noted in the context of time-dependent gravimetry for a purely Newtonian potential [74]. By now considering the cases where the uncertainty in the parameter is a fraction of the parameter itself, we are able to define the regions in which modifications to Newtonian gravity can be established with certainty.

Experimental parameters
We assume that the oscillating source mass oscillates at the resonant frequency of the optomechanical system. We further assume that the source mass is made of solid gold, which has a density of ρ = 19.3 × 10 3 kg m −3 . For a mass of 10 −6 kg (1 mg), this translates into a source mass radius of R S = 2.3 × 10 −4 m. While this mass is very small compared with those currently used in atom interferometry experiments [40], gravitational fields from masses of slightly larger radii have recently been detected [80]. The reason for choosing such a small mass is that the systems can be placed very close together while still achieving a significant oscillation amplitude. This allows us to probe parameter regimes of a short-ranged force. Due to the scaling of the sensitivity as ∆θ ∼ x 2 0 , choosing a smaller x 0 is always going to be beneficial. We therefore set x 0 = 10 −3 m and assume that the oscillation amplitude ratio is = 0.1. This ensures that, when the source mass oscillates, it does not come into contact with the optomechanical system § For the optomechanical probe, we use the following example parameters: we assume that the effective mass of the optomechanical probe is m = 10 −14 kg, and that the light-matter coupling has an amplitude of k 0 /(2π) = 10 Hz. We then assume that the mechanical frequency can be made as low as ω mch /(2π) = 100 Hz, which is important since the expressions for ∆κ and ∆σ scale with ω 5/2 mech . For the squeezed coherent state, we assume that the coherent state parameter is given by |µ c | 2 = 10 6 and that the phase of the squeezed light can be set to ϕ = π, which ensures that the photon number variance (∆N a ) 2 shown in equation (33) is maximized. One of the highest squeezing factors that have been achieved to-date is r sq = 1.73 [109], which is what we choose to include here. We also consider a protocol where we perform M = 10 3 measurements at time t ω mech = 20π, which allows us to improve the sensitivity a bit further.
To derive the bounds on the chameleon parameters M and Λ, we assume that the optomechanical system can be operated in high vacuum. This also helps in terms of mitigating mechanical noise; in generic oscillators, damping effects are well-understood and largely not present below 10 −7 mbar [110]. On the other hand, it can be challenging to confine a levitated optomechanical system at high vacuum [89]. Recently, however, several works have demonstrated trapping at 10 −7 mbar of pressure [89,111], even going as low as 9, 2 × 10 −9 mbar [112]. Using these values as our starting point, we note that 10 −9 mbar translates into a molecular background density of ρ bg = 8.27 × 10 −14 kg m −3 . To derive this value, we have used the ideal gas law, which can be rewritten to give ρ bg = P m N2 /(k B T ). Here, P is the pressure (in Pascal), k B is Boltzmann's constant, T is the temperature (in Kelvin), and where we have assumed that the vacuum chamber has been vented with hydrogen of molecular mass m H = 3.3 × 10 −27 kg before being emptied (that is, it was filled with hydrogen gas, such that any residual particles inside the chamber are H 2 particles).
To see how strong the modified contributions to the force are compared with just the Newtonian part, we plot the amplitude of the time-dependent modification F mod = GmMs . The result can be found in Figure 2, where we have plotted contours for F mod /F N = σ using the experimental parameters in table 1. Figure 2a shows F mod /F N as a function of α and λ, and Figure 2b shows F mod /F N as a function of M and Λ. The filled-in contours in Figure 2b correspond to the force shown in equation (17), where the screening from the optomechanical probe itself has been ignored. The lines, on the other hand, correspond to the force shown in equation (18) where the screening from a spherical probe has been taken into account.

Fundamental bounds for the Yukawa parameters α and λ
We are now ready to compute the bounds on the parameter ranges that could potentially be tested with a quantum optomechanical system. To find the bounds, we consider the ratios ∆κ/κ and ∆σ/σ as functions of α and λ, where κ and σ were defined in equation (24) as the modification due to the gravitational force at the equilibrium distance and the amplitude of the time-dependent contribution. The result can be found in figure 3a: the dark green dashed line shows where the relative error satisfies ∆κ/κ = 1, and the dotted green line shows where ∆κ (mod) /κ = 1. Since κ corresponds to the static modification of the gravitational force, modulating the optomechanical coupling does not improve the sensitivity. We instead focus on the dynamic contribution from σ. The lighter purple region shows where ∆σ/σ < 1, and the darker purple region shows  where ∆σ (mod) < 1. The resonantly modulated optomechanical coupling provides a significant enhancement for ∆σ. The general features in figure 3a can be understood by examining the form of κ and σ, which are shown in equation (24). When λ x 0 , the exponential can be approximated as e −x0/λ ∼ 1. This means that σ becomes σ ∼ 2α, which is independent λ and thereby explains the straight line at |α| ∼ 10 −3 . Once λ < x 0 , which corresponds to a short-ranged Yukawa force, the effect can no longer be seen by the optomechanical probe. However, the bounds in figure 3a could be shifted to the left by decreasing x 0 . Care must be taken that the two systems do not touch, which is limited by the source sphere and probe radii, as well as the oscillation amplitude x 0 . For the example parameters used here, the smallest distance between the system is 0.7 mm.

Fundamental bounds for the chameleon parameters M and Λ
To obtain the bounds on M and Λ, we rescale M in terms of the reduced Planck mass M P . We then compute the bounds for M and Λ by plotting ∆σ/σ as a function of M and Λ for the following two cases: (i) when the probe is approximated as a point-particle (no probe screening), and (ii) when the screening from the probe is taken into account. The latter we denote by ∆σ (scr) and ∆σ (mod) (scr) . We compute these quantities by numerically solving equation (13) for S S and S P at each point. The expression for σ given in equation (28).
The result can be found in figure 3b. Note that we do not plot the bounds for ∆κ and ∆κ (mod) for clarity, and because as static contributions they are more difficult to distinguish from a constant noise floor. The lighter regions show the bounds when the optomechanical probe does not contribute to the screening   [113] and recent results presented in [114] (see figure 6). (b) shows bounds in terms of M and Λ, which are the mass and energy-scale for the chameleon screening mechanism. The experimentally excluded regions are based on those reported in Ref [115].
of the fifth force. This is equivalent to approximating the probe as a point-particle. In contrast, the darker regions show the reduction in sensitivity due to the screening that arises from a spherical optomechanical probe.
To explain the features of the plot, we draw lines where the screening from the probe S P and source system S S vanishes. The magenta line shows where S S = 0 and the orange line shows where S P = 0. Above each line, the screening is zero, while below the lines, the screening lengths increase and the modifications to Newtonian gravity can no longer be detected. Finally, the right-most boundary of the dark purple area can be understood as follows: The appearance of M −3 in m bg (see equation (11)) ensures that, when M is large compared with the other quantities, m bg is small. This, in turn, means that the range of the force λ bg , as shown in equation (15) will be large. It then follows that the amplitude σ (see equation (28)) will be approximately σ ≈ α bg,P , where α bg,P = 2M 2 /M 2 P from equation (19) (note that ξ S = ξ P = 1 because we are considering the range of Λ above the orange and magenta lines). This means that σ is independent of Λ and the boundary becomes a vertical line. The point at which the ratio ∆σ/σ = 1 then occurs is M/M P = 48.1.

Relation to existing experimental bounds
To see how the theoretical bounds in figure 3 relate to known experimental bounds on Newtonian gravity, we plot the convex hull of the shaded areas in figures 3a and 3b against the bounds presented in Refs. [113][114][115].By comparing with experimental results, we are able to demonstrate where optomechanical systems could help further constrain known bounds according to the results in this work. We emphasise however that this comparison is highly hypothetical, since experimental challenges such as noise, long-term stability, and integration over many runs of the experiment have not been included in our analysis. Much more work is required before it is known exactly how the optomechanical probe compares with other platforms (see section 7).
The bounds can be found in figure 4, where figure 4a shows the bounds in terms of α and λ, and where figure 4b shows the bounds in terms of M and Λ. The yellow regions show the convex hull of the bounds derived in this work, and the purple region shows the combined parameter spaces that have been experimentally excluded. The orange area in figure 4b shows the excluded region for when the optomechanical probe is approximated as a point-particle, i.e. the chameleon screening due to the finite size of the probe is neglected.
Our results indicate that, for the values used in this work, even the ideal realisation of a nonlinear optomechanical sensor achieves similar bounds on α and λ to those already reported in the literature. The decoherence, dissipation and thermalisation effects not accounted for in this description are likely to further reduce the sensitivity. This suggests that the sensitivity of the system must be improved further, should we wish to probe the hitherto unexplored regions in figure 4a. From inspecting equations (31), (32), (34), and (35), we note that the strongest dependence is with the mechanical frequency ω m . Thus the lower ω m , the better the sensitivity. Another strategy would be to increase the strength of the light-matter coupling k 0 , however this is a long-standing challenge for many experimental platforms. More effective perhaps would be to decrease the separation distance x 0 between the probe and source systems, which would allow the optomechanical sensor to explore a larger range of λ, in particular smaller λ, since the Yukawa potential will not be as suppressed there. However, as the sensor is moved closer to the source sphere, the Casimir effect is expected to strongly contribute to the resulting acceleration (see below). On the other hand, our results according to figure 4b indicate that optomechanical systems could be used to probe some hitherto unexplored regions of the chameleon parameters M and Λ. The advances here likely depend on the quality of the background vacuum.

Discussion
In this section, we discuss the challenges that must be overcome when considering an experiment of this nature. They include systematics and noise that affect the experiment, as well as forces that arise from the Casimir effect.

Examining the conditions for linearising the force
In order to definitely rule out modifications to the Newtonian potential, we must experimentally determine if the observed data deviates from that predicted by Newtonian gravity. Doing so requires extensive knowledge of the full dynamics of the system, including higher-order contributions from the Newtonian potential that we have neglected in our main analysis. With this in mind, we examine the derivation of the linearised gravitational potential (see the expansion in equation (3)) to determine when this linearisation breaks down. We assumed that the perturbation δx to the position of the optomechanical element is small compared with x S (t) (the distance from the probe to the source mass) at all times. However, depending on the intended precision of the measurement of the force, Newtonian gravitational terms of second order in δx may become relevant, that is, terms of the form ∝ b † +b 2 . These terms can be included into the full dynamical analysis, which has been done in [116]. We leave performing the same analysis for modified gravity to future work. Moreover, the radiation pressure found in an optomechanical setup has the explicit effect of displacing the mechanical element. When the light-matter coupling is modulated at mechanical resonance, the maximum position increases linearly as a function of time [74]. Once this displacement grows too large, the approximation under which the optomechanical Hamiltonian in equation (1) was derived is no longer valid (see e.g. Ref [117] for details of how the optomechanical Hamiltonian is derived). A method for dealing with a displacement driven by radiation pressure would be attempting to cancel the expected radiation pressure by manually introducing a time-dependent linear potential ∼ (b † +b) into the dynamics [74]. In this way, the displacement from the light-radiation pressure is cancelled, while the phase from the gravitational interaction is still imparted on the optical state. The drawback of this method is that it most likely introduces additional noise into the experimental setup from the linear driving term. We do however leave the full quantum metrology analysis to future work.

Limitations due to the Casimir effect
Due to the relative weakness of gravity compared with electromagnetic force, the latter are likely to dominate any experimental setting. Therefore, any stray electromagnetic effects must be controlled very precisely in order to detect deviations from Newtonian gravity. One of the most important effects that has to be taken into account is the Casimir force [118], which becomes significant when the distance between the probe and the source mass is small. To estimate the effect of the Casimir force, we use an analytic formula given in [119] (based on the results of [120]) for the force due to the Casimir effect between two homogeneous perfectly conducting spheres at a distance much larger than their radii. The model of two perfectly conducting spheres is unlikely to accurately describe the experimental realisation of optomechanical setup described in this article, both in terms of geometry and material. Therefore, we will use this case to give only a first estimate of the effect and discuss how to suppress it.
We consider the Drude boundary condition model for isolated conductors (see [119] for details). For the distance between the probe and the source x 0 − R S − R P being much larger than the thermal wavelength, i.e. x 0 − R S − R P λ T = c/(2πk B T ) (where the thermal wavelength is about 1µm at room temperature) the classical thermal contribution to the Casimir force dominates, which leads to the expressions where R P and R S are the radii of the probe and the source, respectively. At room temperature and for the parameters given in table 1, equation (36) leads to an acceleration of the order of 9 × 10 −13 m s −2 , experienced by the probe mass, while the gravitational acceleration induced by the source mass is of the order of 6×10 −11 m s −2 . Casimir forces are therefore of order 10 −2 smaller than the main gravitational component. The size of the fifth force corrections we consider here are largely controlled by the two parameters σ and κ as shown in equation (25). At peak sensitivity, when x 0 ∼ λ, this means that the ratio of the Casimir force to the Newtonian gravitational force should be compared to α. We see from figure 3a that this ratio can be as low as 10 −3 at the edge of the detectable region; the corrections are even smaller at lower λ. Furthermore, since the Casimir force grows very strongly with the inverse distance of the source and probe mass, the Casimir force quickly overshadows the fifth force contributions by many orders of magnitude when the source-probe distance is decreased to achieve better sensitivities. This shows that the Casimir effect is a relevant systematic that has to be controlled, that is, either precisely quantified or reduced. One way to reduce the force is to lower the temperature of the setup. Another option to suppress the Casimir effect is to place a material in-between the source mass and the sensor that acts as a shield to the Casimir effect [121,122]. The Casimir force of the shield will be stationary while the un-shielded gravitational acceleration will be time-dependent, and therefore, clearly distinguishable [123]. This approach is, however, limited by the size of the shield. For example, in levitated optomechanics, the screening scheme can be naturally realized by placing the source mass behind one of the cavity end mirrors such that the mirror serves as a shield. However, in the case of detecting modifications due to a chameleon field, the presence of the mirror might introduce additional screening effects that need to be accounted for. The Casimir effect may also be reduced by modulating or compensating for the Casimir force with radiation pressure [124], nano-structuring of the source and probe surfaces [125], or an optical modulation of the charge density [126].
Further analysis of the impact of a shield, or other techniques for accounting for the impact of the Casimir force, will require detailed numerical modelling. For example, Pernot-Borr et al. [78] considered the impact of cylindrical walls on the screening of a source, finding that it can depend strongly on the thickness of the wall used for screening. Since we here consider the fundamental limits of an optomechanical setup, we leave a numerical analysis of the impact of different approaches to future work.

Improvements to the sensitivity
There are a number of ways in which the sensitivity of the optomechanical system can be further improved. In this work, we considered spherical source masses and probes in order to analytically derive the screening from the probe, however, choosing a different shaped source may improve the bounds that could be achieved. For example, a source mass in the form of a slab much larger than the probe system would mitigate gradient contributions from the Newtonian part of the potential, since the gravitational force from an infinite plane is constant. Furthermore, it was shown in Ref [100] that symmetric source masses tend to be much more strongly screened (and thus have smaller detectable effects) then asymmetric sources. Therefore, we would expect to obtain more favourable precision bounds than those presented in this work by considering asymmetric sources. An interesting prospect also arises from the fact that the optomechanical probe itself can also be asymmetric, e.g. in the shape of a levitated rod [127], which offers an additional avenue compared with, for example, atomic systems. However, these non-spherical cases bring with them additional challenges. The approximation used in equation (12) assumes a spherical source (probe), and approximates the nonlinear solution of the chameleon field equation with an analytic expression derived by asymptotic matching. To accurately obtain measurements with a non-spherical setup would require precise numerical modelling of the chameleon field around a (non-spherical) source and probe, such as done in Ref [100]. The precise effect on the sensitivity is left to future work.
As a final note, we mention that the nonlinear radiation-pressure term in the Hamiltonian in equation (1) appears in many different contexts, of which not all fall under the category of optomechanics (such as for example electromechanical setups [128]). Our results therefore apply to these systems as well. We therefore have a large range of systems to choose from when it comes to optimising the geometry and resulting sensitivity.

Future work towards an experimental proposal
The sensitivities calculated in this work give us an indication of the resolution of the force that the optomechancial probe can achieve in principle. That is, we learn the magnitude of the gravitational force that can be detected. In practice, however, we must then determine whether this force is simply the Newtonian force, or whether it is due to the Newtonian force and an additional force that arises from the modification. With a good-enough resolution, such a modification can be detected even if the Newtonian force is much stronger than the modification. There are several methods by which the modification can be detected. The first is to very carefully model the influence of the Newtonian force on the optomechanical dynamics and data that is collected through e.g. a homodyne measurement. If a deviation in the collected data is then seen, steps should be taken to rule out any other source. Another way is to carefully change the equilibrium separation distance x 0 between the source sphere and the optomechanical probe. Since the modifications considered in this work changes quite drastically due to the inclusion of the exponential term in equation (7), it should be possible to detect such a exponential change in the data. Both of these methods here can be theoretically explored in future work.
Our results can be used to evaluate the fundamental ability of a quantum optomechanical system to probe a particular parameter regime of modified gravity theories. A realistic optomechanical system, however, will be affected by a number of systematics and noise sources, including optical dissipation from photons leaking from the cavity, mechanical thermal noise, Brownian motion noise, damping effects, and noise from the trapping or clamping mechanism, as well as radiation back-action noise and shot-noise. Yet additional noise sources include external gravitational noise and environmental vibrations (see e.g. [81,123] for a discussion of a related experimental setup). Generally, such noise sources have spectral contributions at the resonant frequency of the sensor and are enhanced as well as the signal from the source mass that we wish to detect. Therefore, in practice, it may be favourable to consider an off-resonant sensing scheme, such as those discussed in Refs [81,123]. We also note that such additional noise sources will be particularly dominant when the mechanical frequency is low, however we see from equation (31) and (32) that a low mechanical frequency is a necessary requirement if we wish to achieve a high sensitivity. We also note that it is not clear how the sensitivity gained from e.g. modulating the optomechanical coupling changes when the Q-factors of the cavity and the oscillators are considered.
To model the noise and systematics mentioned in the previous paragraph, a plausible next step beyond this work involves linearising the optomechanical dynamics around a strong coherent input-state [59]. With the help of phase-space methods [129], it is then possible to include most of the systematics and noise terms mentioned above into the dynamics. In addition, a homodyne measurement could be modelled using inputoutput theory for the optical mode. One can then examine the susceptibility of the mode and determine the noise levels required for these effects to be detectable [130]. An important question that must be addressed is the laser power required to maximise the sensitivity. Since the linearisation gives rise to equations of motion that differ from those used here, it is difficult to predict what the resulting bounds on modified gravity theories will look like compared with those presented here. Most likely, the presence of noise and absence of non-Gaussian resources (which arise from the nonlinear coupling) means that the prediction for the sensitivity is reduced.
To instead extend the analysis in this work even further in the nonlinear optomechanical regime, we must include noise in the solution of the dynamics for the nonlinear Hamiltonian in equation (1). However, since the resulting nonlinear Langevin equations are generally much more difficult to solve (although certain solutions in the weak-coupling limit and for systems with weak optical decoherence exist [131,132]), we expect this to be challenging. A preliminary step towards modeling Markovian optical decoherence affecting the intra-cavity state was recently taken [133], and mechanical thermal noise has been modelled using a range of methods [134,135]. For a strongly coupled system, however, optical and mechanical noise cannot be treated separately, and must instead be considered together [136,137]. To our knowledge, fundamental quantum metrology bounds in the noisy nonlinear regime have not yet been considered.
Another aspect that needs to be modelled is the additional screening that arises from the inclusion of a shield to block out Casimir forces. In addition, for a levitated optomechanical sphere, a mirror must be placed between the optomechanical probe and the source, which also contributes to the screening (but which can, at the same time, act as the Casimir shield). To carry out a full analysis of the screening, the geometry of the vacuum chamber, along with the trapping mechanism of the optomechanical system and the Casimir shield, must be carefully modelled. It is then possible to exactly predict the magnitude of the modified force that the optomechanical probe can detect.

Conclusions
In this work, we derived the best-possible bounds for detecting modified gravity with a quantum optomechanical sensor. We modelled the effects of a force from an oscillating source mass on the optomechanical probe and estimated the sensitivity of the system by computing the quantum Fisher information. In particular, we considered the additional screening that arises due to the relatively large size of the optomechanical probe. Our results show that optomechanical sensors could, in principle, be used to improve on existing experimental bounds for the chameleon screening mechanism, although more work is needed to evaluate the prospects for using experimental optomechanical systems as probes for modified gravity.

Data availability statement
The code used to compute the screening and sensitivity to chameleon fields can be found in the following online GitHub repository: https://github.com/sqvarfort/modified-gravity-optomech.

Appendix A. The chameleon mechanism
In this appendix we briefly review the derivation of the chameleon mechanism and how it gives rise to a fifth-force; the reader is directed to Refs [12][13][14] for further details. Throughout this appendix we will use energy units ( = c = 1) for notational simplicity. The basic idea of the chameleon screening mechanism is to screen the effects of additional degrees of freedom in a modified gravity model (typically light scalar fields), by making their mass dependent on the local density. This results in a scalar field whose mass is large inside the solar system where the average density is high and is thus difficult to create in collider experiments, but has a lighter mass in the intergalactic medium where the density of matter is lower. Typically, this is achieved using a scalar field whose action is of the form: where g µν is the spacetime curvature, R is the Ricci tensor, V (φ) is the chameleon potential, and S (m) is the matter action. Various choices of screening mechanism are possible with this action, but the Chameleon mechanism corresponds the following choice of the functions V (φ) (the interaction potential) and Ω(φ) (which represents direct, non-minimal coupling between the Chameleon field and gravity): In this work, we will consider the case n = 1 for simplicity. One way to understand the effect of the Ω(φ) term is to regard matter as coupling to the so called Jordan frame metric,g µν = Ω −2 (φ)g µν , while the Chameleon field sees a different metric, g µν , which suggests quanta of the scalar field, if we were somehow able to isolate them, would be observed to fall differently to normal matter (violating the equivalence principle). One can either regardg µν as the "real" metric, in which case φ has an unusual direct coupling to gravity, or g µν , in which case all particles have a special coupling to φ via the function Ω(φ) appearing wherever the metric does in the matter action, S m . Ultimately what matters, however, is how objects will be observed to move in the presence of this scalar field. Formally, we can obtain this from the geodesic equation for the 4-vector position X µ = (t, X), since normal matter sees the metricg µν : If we regard g µν as the true metric, then the effect of the chameleon field is to add what appears to be a fifth force, since when we take the Newtonian limit we can re-writeΓ ρ µν in equation(A.4) in terms of the Newtonian potential, Φ N and Ω to obtain: where Φ C = − log Ω(φ(X)) is an effective fifth-force potential. The strength of this fifth force is characterised by Ω, but to compute its effects we need to know how the scalar field couples to matter. Varying equation (A.1) with respect to φ we obtain: where T µν is the Hilbert Stress energy tensor. For non-relativistic matter, the scalar field is sourced by the local matter density, with g µν T µν = −ρ, giving This is equivalent to the scalar field interacting via the potential: If we use the form equation (A.3), then for cases where φ M , we can approximate log Ω(φ) as: Under this approximation, the effective scalar field potential becomes (A. 10) and the Chameleon fifth force potential is: For the purposes of this work, we will consider the n = 1 chameleon field. In a region with constant mass density ρ bg , this means that the Chameleon rests at the vacuum value: for which fluctuations of the field have mass: As expected, the mass of field fluctuations increases with the background density, which means in areas of comparative high density such as inside the solar system , the mass is large and the scalar field difficult to excite and detect.

Appendix A.1. Chameleon Field From a Spherical Source
The chameleon field in the vicinity of a spherical source of mass M S and radius R S can be computed by solving the Klein-Gordon equation, This is a non-linear equation, but an approximate solution was found by Burrage et al. [33] in the limit m bg R S 0. This is valid for the atoms considered there, but in our case we may need to consider larger sources. We therefore repeat the derivation of Burrage et al. [33] for the case of arbitrary m bg R S . The fundamental strategy uses the method of asymptotic matching [138] to derive an approximate solution over for the full domain of the differential equation, by smoothly matching together solutions valid in different domains.
In the strongly perturbing case (ρ S R 2 S > 3M φ bg /( c)), the solution reaches its equilibrium value, φ S , for density ρ S at some radius S ≤ R S . We denote the region with r < S the interior region, or region I. The solution there can thus be approximated as constant There is then a transition layer (region II) between S < r < R S where the solution rapidly shifts towards the background value, φ bg . Since the density ratio between the source-sphere and the external vacuum is high, we will find φ bg φ S as a result of equation (A.12). In the transition layer (S < r < R S ), the field Compared to the average density inside a cosmological void.
will begin to increase, eventually reaching a regime where φ φ S . Because we can re-write equation (A.10) as 16) then once φ φ S the density-dependent term dominates the potential. Under such conditions, V eff (φ) ≈ ρ/M and we can solve equation (A.14) analytically: Finally, far away from the source sphere, φ is close to its background value, φ bg , and we can approximate the potential as quadratic: Here, region III is defined as r > R S . Note that although φ φ S outside the sphere, equation (A.17) does not apply because the density outside the sphere is now ρ bg which is typically much smaller than ρ S : the density dependent term in the potential is thus no longer dominant outside the source. We note, however, that solutions φ II and φ III are technically only valid in the vicinity of r ∼ R S and r R S respectively. However, we can approximate the behaviour of the fully-non-linear solution for all r by matching these asymptotic solutions at S and R S , which imposes four constraints to ensure smoothness of the asymptotically matched solution: We also note that we require F = 0 to have a solution approaching φ bg as r → ∞, which means that there are four unknowns, C, D, E, and the radius S. We solve for these four unknowns, finding Taken together, Eqs (A.19)-(A.21) and (A.22) imply equation (12), and reduce to the Burrage et al. [33] result in the case m bg R S → 0 limit. As a cubic equation for S, it is of limited use to express S in closed form, and for the purposes of this work we solve equation (A.22) numerically. This can run into catastrophic cancellation problems, due to the finite floating-point precision, when S ≈ R S (that is, the heavily screened regime), since we need to compute 1 − S 3 /R 3 S . Hence, when the numerical solution gives S/R S close to 1, we switch over to an analytic approximation obtained by substituting S/R S = 1 + and solving for to first order. This gives We see that equation (A.24) agrees with the Taylor expansion of equation (14) in the m bg R S → 0 limit.

Appendix B. Time dependence of a chameleon field
The main result of this work is that an oscillating optomechnanical system can be used to detect the presence of chameleon fields. However, throughout we have made the assumption that such a field responds essentially instantaneously to the motion of the mass that sources it. For purely gravitational effects this assumption is justified since information about the position of the source mass propagates outwards at the speed of light. However, the situation is less clear for a scalar field sourced by a moving mass, since the field is massive: there is no a priori guarantee that information can propagate outwards at the speed of light. To address this, we derive the propagator for a chameleon field sourced by a point mass and demonstrate that the resulting retarded chameleon potential can be treated as if information propagates instantaneously, for non-relativistic oscillating masses.
Appendix B.1. Time dependence of the gravitational potential First, we will derive the response of the gravitational field to a small mass, neglecting back-reaction and gravitational wave emission as negligible. In the linearised limit which allows us to make contact with Newtonian gravity, the metric perturbation h µν around η µν satisfies: where we are using the − + ++ metric convention. We choose the Lorenz gauge, defined by the condition ∂ µ h µν = 1 2 ∂ ν h, which simplifies this to: whereh µν = h µν − 1 2 hη µν is the trace-reversed perturbation. A generic expression for the Hilbert stress energy tensor is: where S m is the matter action. For a point particle of mass M S , the appropriate matter action is: where q µ (τ ) describes the particle trajectory and dots denote differentiation with respect to proper time, τ (the variation with respect to q µ yields the geodesic equation, verifying that this is indeed the action we seek). This means that the stress-energy tensor at position X for a point particle following trajectory q(τ ) is: For low (non-relativistic) velocities in a Minkowski background this means T 00 = M S δ (3) (X − q(τ )) = ρ, as we would expect, with all other components zero. Note that we are ignoring any back-reaction from this moving particle, which would give only higher order corrections. Taking the trace of equation (B.2) yields: which is the same equation satisfied byh 00 , and so we have h =h 00 (note the sign -η µνh µν = −h, but η µν T µν = −ρ). Meanwhile, ∂ λ ∂ λh ij = 0 as there are no spatial parts of T µν . This describes propagating gravitational waves, which we neglect, and can be safely set to zero. Thus, we need only solve equation (B.6), which has a known solution in terms of a retarded potential: For the point source with density ρ(X, t) = M S δ (3) (X − q(t)), this means: where t ret is the retarded time that solves: Note, comparing with the standard perturbative parameterisation of the metric: we see that Φ = −h/4, or in other words: which is the retarded gravitational potential for a point source of mass M , as we would expect. Note that this is well defined (in the Lorenz gauge), as it is not derived from energy considerations but from the equation of motion of the metric perturbation.
We use the formula: where X i are solutions of f (X i ) = 0. In the q constant case, it is easy to see that there are two solutions, t = t±|q−X|, but only the negative solution matters, due to the θ(t−t ) term (this is the causal effect of the retarded propagator, and ensures that we only integrate over contributions from the past of the time t we are looking at). In the case where q is time dependent, solving (t − t ) 2 − |X − q(t )| 2 = 0 is less trivial, but still results in a unique retarded time, t ret , exactly the same quantity that is well known from electrodynamics. To see that it is unique, consider that t ret is by definition the time at which light arriving at an observer at time t was emitted by the source. Assume there are two such times, t 1 , t 2 . We have (temporarily putting back the factors of c for clarity) c(t − t i ) = |q(t i ) − X|, and can subtract these two equations from each other to obtain c(t 2 − t 1 ) = |q(t 1 ) − X| − |q(t 2 ) − X|. This implies that the distance of q from X has changed at the speed of light -not possible unless the source itself is moving at the speed of light. Thus, for sub-luminal sources, t ret is unique ¶. Hence: where we have used equation (B.9) and v(t) is the velocity of the source. Thus, the first part of the integral reduces to: .

(B.20)
This is indeed exactly the same as the expression found in electrodynamics, where the propagating field (the photon) is massless, and thus the second term in equation (B.17) is not present. In our case, however, we have to deal with the massive part of the integral too: In this case, the effect of the two Heaviside step functions is to force us to integrate over the past, up to the retarded time: Now, make the substitution: to obtain: This would not seem to offer a significant simplification, unless we adopt the low-velocity approximation, |v 0 |/c 1. The second term in the denominator is O(v/c), so we can in general neglect it (note that ¶ Actually, there is still a solution for t adv > t, but this is eliminated due to the causal θ(t − t ) function.
|X − q(t )|/c is not in general small as X can be arbitrarily far from the source). We can also expand |X − q(t (u))| 2 as a power series in u around the retarded time, t ret = t (0): We find: , (B.28) And evaluated at u = 0 (or t = t ret ), this gives: .

(B.33)
Hence: We substitute in equation (B.9) to obtain The u 2 term is proportional to v/c once SI units are restored and so we can ignore it in the non-relativistic limit. Generally speaking, higher order terms in the expansion about u = 0 will also have terms proportional to v/c 1, so we neglect them + . This reduces the integral to: (B.34) + One can argue that this must be true for the solution to reduce to the static Yukawa potential in the v → c limit.

(B.35)
This is the expected Yukawa potential, only with the retarded time for the position of the source. It is worth noting that the fact that the force carrier (in this case the scalar bosonic chameleon field excitations) is massive does not appear to affect the retarded time, which describes information propagating through the field at the speed of light, even though the bosons themselves do not. The fact that the retarded time appears, both here and in electrodynamics, is not because the force carriers themselves (photons in the case of electrodynamics) travel at the speed of light -it is due to the causal structure of space-time itself.

Appendix B.3. Screening by the optmomechanical Probe
Burrage et al. [33] derived an expression for the force between two extended spheres due to a chameleon field. However, their approach considered the forces between individual atoms, for which the range of the force could largely be ignored (corresponding to the m bg r 1 limit). With the larger sensor devices we consider in this work this is not necessarily applicable. Furthermore, since in this case we envision an oscillating source, we need to be certain that time-dependent effects do not come into play. For this reason, we re-derive the force between two spheres without the static spheres and m bg r 1 assumption. The derivation closely follows that of Burrage et al. [33], with the assumption of a spherical probe to provide a simplification (more accurate modelling of the optomechanical probe will in principle be necessary for performing the actual experiment).
As in Burrage et al. [33], we denote the source sphere A and the test sphere (B) for which we are computing the force (as a model for the optomechanical probe). The force is determined by the rate of change of the momentum resulting from the energy momentum flux across the surface of ball B: where τ j i are the spatial components of the energy momentum tensor (including gravity), ∂B is the surface of ball B, n j the surface normal vector and dS the surface area element. The gravitational and matter contributions to the energy momentum tensor follow as in Burrage et al. [33], with the main change being to the chameleon field term: (B.37) The chameleon field at position X and time t with X centred on ball B is given by (B.38) where r = |X| and where the prefactors ξ A and ξ B are given by Note that that we have assumed the linear regime, where the fields from the two chameleon sources add. The gradient of the contribution φ B (x) from ball B is given by The chameleon stress-energy tensor is therefore given by whereas Burrage et al. [33], we ignore the contribution from the potential, V (φ). Note that Burrage et al. [33] assume the m bg r 1 limit, which we generalise here to arbitrarily large r. To do this, we will make only the assumption that |X|/|X A | 1 on the surface of ball B, that is, that the radius of ball B is much smaller than the distance to ball A. The means we can expand We can then expand the derivative of the field from sphere A as where we have deliberately not expanded the exponential in |X|/|X A |. This is because we cannot generally assume that m bg |X| << 1 (though we will examine this limit later). We are seeking to perform an integral of the form First, use the fact that ∂ t φ B = 0 (we are in the frame of reference of B, so all time dependence is in the motion of ball A). The time derivative of the field is which can be written in terms of the spatial derivative because the time dependence only appears through X A , which always appears together with X in the form |X − X A |. This allows us to simplify the expression for the stress-energy tensor: We therefore find that the integral splits into three parts: (B.52) Similarly: This means we can write the three integrals I AA , I AB , I BB as (B.57) where we have expressed everything in terms of the following integrals: In fact, one can immediately see that I BB = 0, because the integral of a cubic polynomial vanishes over a sphere: This follows from symmetry since the integrand is odd when X → −X and we integrate over the whole sphere (one can also verify this directly by evaluating it in spherical polar co-ordinates).
To evaluate the remaining integrals in equations (B.59)-(B.61), we switch to polar co-ordinates. Without loss of generality, we can arrange our co-ordinates at any time to be such that X A is along the z axis, which simplifies the calculation. For B i (m) we obtain We therefore only need to evaluate B 3 : We can summarise this as which generalises, for an arbitrary axis choice, to Next, we evaluate B ij (m). We can make similar arguments here to conclude that B 13 = B 23 = 0 (for the same reason that B 1 = B 2 = 0). Furthermore: Hence, all the cross-terms vanish, so this tensor is diagonal. The diagonals are: Finally, we can compute B ijk (m). First, note that this is a rank-3 symmetric tensor. Rank m symmetric tensors in d dimensions have (m + d − 1)!/(d − 1)!m! independent components, which for d = 3 gives (m + 1)(m + 2)/2. This means we have to compute 10 independent components in total. However, as before most are actually 0. Again, similar arguments to before imply B 331 = B 221 = B 321 = 0. We also find that: To summarise this, we note that only terms with at least one of i, j, k equal to 3 are non-zero. The other two indices must then be a diagonal matrix. This means that we can write: +δ j3 δ k3 πR B m 4 (4mR B (9 + m 2 R 2 B ) cosh(mR B ) − 4(9 + 4m 2 R 2 B ) sinh(mR B )) + (all permutations). (B.79) Or, using symmetrised index notation: We can now proceed to substitute these into the expressions for I AA and I AB . First, we consider I AA , since this term is argued to be zero by Burrage et al. [33] in the m bg r 0 limit. First, we can ignore the velocity dependent terms, since we work in the non-relativistic limit where |v A | 0 (we are using units where c = 1). To gain some further understanding of the behaviour, let us consider the R B → 0 limit: This gives

(B.88)
This suggests that I AA is subdominant in this limit, provided: Then this means that we can essentially ignore the I AA term in the R B → 0 limit, which agrees with the calculation of Burrage et al. [33]. Furthermore, provided we stay in the limit in equation (B.89), then I AA will always be suppressed relative to I AB . Generally speaking, this may not be the case, however, in which case we would have to include the I AA term. For now, let us compute the force without this term, as just include the I AB contribution: In the same co-ordinates, we find which implies: (B.95) We conclude that (m bg R B cosh(m bg R B ) − sinh(m bg R B ) .

(B.96)
Including the gravitational force, and moving to co-ordinates centred on ball A, we find: where we have a form-factor function that modifies the force, given by This force has several peculiar features. Firstly, there appears to be a double-counting of the distance from the centre of ball B to its surface (R B ) featuring in the exponential suppression with distance. Note, however, that this additional exponential arises from the screening affect of the probe itself, that is, its lack of response to the chameleon field when the probe is itself large. Secondly, we note that if we computed the force of ball B on ball A, we would not obtain a symmetric force of opposite sign. In other words, this force appears to violate Newton's third law, which on the face of it suggest that momentum is not conserved. However, this is misleading, because the derivation takes into account the stress energy tensor not only of the ball itself (τ (m)µν ) but of the chameleon field (τ (φ)µν ). Momentum should still be conserved if the momentum of both the balls and the field are included in the calculation. This differs from the Newtonian force derived by Burrage et al. [33] since in the m bg R B 1 limit the force is inverse-square, meaning that Newton's shell theorem applies; a sphere should exert the same inverse-square force as an equivalent point mass. But this does not apply to Yukawa potentials, implying that the chameleon force does not obey the where the coefficients are given by For an initially coherent state in the optical field and a thermal state of the mechanical element, as that shown in equation (2), we have λ n = tanh 2n (r T )/ cosh 2 (r T ) and |λ n = |ζ ⊗ |n . The QFI can then be written as the following expression given the initially coherent state of the optical mode and thermal state of the mechanical element shown in equation (2) I θ = 4 B 2 |µ c | 2 + 4 We can then obtain the expressions in equations (31), (32) and (34), (35) by dropping the second term, which we can assume to be much smaller than the first term, and taking the derivative of FN a and FB + for κ and σ while assuming k(t) = k 0 and k(t) = k 0 cos(ω m t), respectively.
We emphasise here that the reason that the time dependent Newtonian gravitational acceleration does not appear in the result is due to the linearity of the derivative. Ultimately, a sensing scheme of this form has a certain resolution, which we are able to compute from these results. In practise, the data must still be analysed in order to distinguish between the Newtonian and the modified gravitational force, which we discuss in section 7.