Dynamical response of Bose-Einstein condensates to oscillating gravitational fields

A description of the dynamical response of uniformly trapped Bose-Einstein condensates (BECs) to oscillating external gravitational fields is developed, with the inclusion of damping. Two different effects that can lead to the creation of phonons in the BEC are identified; direct driving and parametric driving. Additionally, the oscillating gravitational field couples phonon modes, which can lead to the transition of excitations between modes. The special case of the gravitational field of a small, oscillating sphere located closely to the BEC is considered. It is shown that measurement of the effects may be possible for oscillating source masses down to the milligram scale, with a signal to noise ratio of the order of 10. To this end, noise terms and variations of experimental parameters are discussed and generic experimental parameters are given for specific atom species. The results of this article suggest the utility of BECs as sensors for the gravitational field of very small oscillating objects which may help to pave the way towards gravity experiments with masses in the quantum regime.


Introduction
Bose-Einstein condensate (BECs) are very small and extremely cold systems of a large number of atoms. These properties are famously exploited for high precision measurements of forces using atom interferometry [1][2][3][4][5][6].
Another method that utilizes BECs as sensors for forces is to measure the forces' effect on the collective oscillations of the atoms in the BEC. One specific example is the measurement of the thermal Casimir-Polder force [7], which was theoretically proposed in [8]. Further theoretical proposals to use collective oscillations in BECs and, in particular, their phonon modes for sensing purposes include, for example, gravitational wave detectors [9][10][11] , sensors for the effect of spacetime curvature on entanglement [12] (for a review see [13]) and magnetic field and rotation sensors employing solitons formed by optical lattices [14].
In this article, we investigate the effect of an oscillating gravitational field on the phonon modes of a BEC in a uniform trapping potential for the particular case of the gravitational near field of a small, oscillating gold or tungsten source mass. In particular, we show that the effect may be detectable for source masses down to the milligram scale considering state of the art experimental parameters. The abilities of experimentalists to cool and control BECs of large numbers of atoms are advancing quickly and BECs may become sensors for very weak, oscillating gravitational fields in the future. The particular situation of a small source mass could be used to measure the gravitational field of a macroscopic quantum system in a spatial superposition state. The preparation of such a state is proposed in [15] and a proof of principle experiment for the measurement of the gravitational field of a small, oscillating mass was proposed in [16], which is based on a macroscopic test mass instead of a BEC. Furthermore, the experimental situation considered in this article may be useful to approach the experimental realization of the proposal to use the phonon modes of a BEC to detect gravitational waves presented in [9] and the proposal to measure gravitationally induced Josephson tunneling presented in [17]. While a relativistic framework is used in [9,17], we consider only the effect of a Newtonian potential in this article. However, in [18], it was shown that the effect of a gravitational wave on the detector proposed in [9] can be completely described in the Newtonian limit. Additionally, a Newtonian potential can be derived from a spacetime metric using the proper detector frame [19] in a similar way to the investigations in [20] for an optical resonator in curved spacetime. Furthermore, the equations describing the BEC in this article can be derived from the non-relativistic limit of the relativistic equations used for the derivations of the results in [9,17]. Therefore, the results presented in this article can be used to describe the measurement of gravitational effects with BECs in a general relativistic context using the techniques employed in [20].
Besides the above mentioned ambitious benefits of an experimental realization of the situation presented in this article, such a realization would be the first investigation of the interaction of an oscillating gravitational field with the phonon modes of a BEC, which may be considered an achievement by itself.
We start our considerations by discussing our approximate description of the gravitational potential seen by the BEC in section 2. In section 3, we introduce our first description of the BEC in terms of the order parameter and we introduce the basic equations for perturbations of the order parameter which represent the phonon modes. Then, we introduce dissipation in the phonon modes in section 4, and we connect the amplitude of the perturbation of the order parameter to the number of created phonons in section 5. In section 6, we investigate the effect of resonant driving by the gravitational near field of the source mass on the BEC. In section 7, the quantum field theoretical description of phonons is introduced and we derive the amplitudes for particle creation due to the oscillating sphere. In section 8, we derive expressions for experimental parameters required to achieve a given signal to noise ratio. Furthermore, we give a list of generic experimental parameters that may allow a detection of the phonons created by the gravitational acceleration and the gravity gradient, respectively, due to an oscillating massive sphere of mass M for the two cases of M=200 g and M=0.2 g. In section 9, we conclude and discuss the possibility to use phonons in a BEC to detect oscillating gravitational fields and the influence of seismic and Newtonian noise on the noise background.

The gravitational field of an oscillating sphere
We assume that the BEC is of length L in the oscillation direction of the sphere (see figure 1) and the trap potential is a uniform box, as was realized in experiments such as the one described in [21]. Furthermore, we assume that the BEC is much smaller than the source mass, so that we can restrict our considerations to one spatial dimension. Let us denote the distance from the center of the trap potential in the direction of the massive sphere as x and let R(t) be the distance between the center of the trap potential and the center of the sphere. For small L, the Newtonian potential of the moving sphere can be approximated, for all positions x inside the trap potential, as 3 are, respectively, the value of the gravitational potential at the center of the box trap potential, the gravitational acceleration due to the sphere evaluated at the center of the trap potential 3 and its first derivative in the x-direction, the gravity gradient. By comparing the second and third term in equation (1) in the range of the BEC, we find that the second term is always larger than the third term by a factor ( ) R t L 2 . Therefore, ( ) G t becomes significant only for distances from the center of the sphere that are a few orders larger than the extension of the BEC trapping potential. When we investigate the perturbations of the BEC induced by the Newtonian potential (1), we will find that the gravity gradient leads to a very different signature than the acceleration. This is because the second term in equation (1) is linear in x, while the third term is quadratic in x.
Let us assume that the source mass moves sinusoidally about a fixed position so that . If we assume that R 0 ?δ R , we can expand a(t) and ( ) G t in δ R . In the following, we will assume that δ R /R 0 is small enough to stop the expansion after the linear order in δR. We find  In an experiment it may be possible to move the mass such that the acceleration or the gravity gradient become exactly sinusoidally time-dependent. Hence, the assumption of small values of δ R /R 0 is a restriction that may be lifted by a different ansatz and more detailed calculations. For m = L 200 m and δ R =2 mm, a minimal distance between the end of the box potential and the surface of the source mass R min =1 mm and a 200 g tungsten or gold sphere, we have R 0 =r+R min +L/2+δ R ≈17 mm, where r≈14 mm is the radius of the sphere. We find that the amplitudes of the terms proportional to In the next section, we will introduce the description of the BEC and show how to derive the effect of the oscillating terms in the gravitational potential.

BEC mean field perturbations
One way to describe a BEC and the perturbations induced by an external potential is in terms of the order parameter ψ governed by the GP equation and perturbations governed by the Bogoliubov-DeGennes (BDG) equations derived from the GP equation. A second option is a microscopic description in terms of the atom field which is split into the condensate fraction and a phonon field. The condensate fraction is again described by the GP equation. The phonon field evolution can be given by an interaction Hamiltonian and a mode decomposition with modes governed by the BDG equations. In sections 3-5, we follow the first option. The second option is followed in section 7. Both approaches are only valid for BECs at low temperature. In particular, we need m  k T B , where μ is the chemical potential of the condensate, T its temperature (see section 8 for explicit values of these parameters) and k B the Boltzmann constant. For , only a small part of the atoms is not in the ground state and its back action on the order parameter can be neglected for the purposes of this article 4 . In this article, we include effects of finite temperature by using an effective description in terms of dissipation in the time evolution of phonons in section 4. Furthermore, we consider effects of evaporating atoms in section 8.
Since the gravitational potential (1) only depends on x, we only consider modulations of the order parameter in the x-direction and use an effectively one-dimensional description. Then, the GP equation is given as [23]    y l y y where m is the mass of the BEC atoms, λ=8π a scatt describes the interaction between the atoms, a scatt is the scattering length and V=V trap +mΦ is the external potential consisting of the trap potential V trap and the Newtonian potential (1). To gain some intuition for the meaning of ψ, we can use the Madelung representation of the order parameter y r = q e i , where ρ can be identified as the atom number density of the condensate and ρÿ∂ x θ/m gives the probability current in the x-direction (see equation (5.10) of [23]), which can be interpreted as the spatial flow of atoms, while ∂ t θ is proportional to the chemical potential μ of the condensate in the case of a stationary solution.
The effect of the gravitational potential of the source mass can be expected to be small. Therefore, it will lead to small perturbations of the atomic cloud. The equations that govern these perturbations can be derived from the GP equation [23,24]. To this end, we set y r d y = + q f , where we assume that r q e 0 i 0 is a solution of the GP equation (4) for the case where the gravitational potential vanishes and the density is stationary ∂ t ρ 0 =0. The complex function δψ is a space and time-dependent perturbation that corresponds to phonons and f 0 is a spatially constant, time-dependent function that captures the perturbation of the ground state energy. We rewrite the GP equation (4) for the unperturbed solution r q e 0 i 0 as the pair of coupled partial differential equations 4 The effect of the thermal cloud of atoms can be taken into account by a generalized GP equation using the Popov approximation. For more details, see for example [22] and section 13 of [23]. In particular, the back action of the thermal cloud leads to a stationary deformation of the order parameter, which would result in a perturbation of the number of created phonons. Since the number of atoms in the thermal cloud is very small in comparison to the number of atoms in the condensate fraction, the effect will be only a small change of the number of created phonons which we will neglect.    q r r q lr r r r q r q where the prime and the dot denote the derivative in the x-direction and the time derivative, respectively. We consider a uniform box potential with infinite potential walls, which implies the boundary condition r = 0 0 at the potential walls. We assume that the length of the condensate L is much larger then its healing length z l r = 1 0 . Then, the density of the BEC can be assumed to be constant everywhere except for a small region close to the potential walls [21,23,25]. We can assume that ρ 0 is constant over the whole length of the uniform trap by additionally restricting our considerations to perturbations δψ with wavelength much larger than the healing length. In particular, this means that we do not consider quasi-particles that behave like free particles.
Inside the box potential, it follows from (5) and r ¢ = 0 0 that q = 0 0 . At the boundary ρ 0 vanishes, r ¢ 0 becomes singular and q ¢ 0 has to vanish for the second equation in (5) to be fulfilled. Since q = 0 0 inside the box potential, it follows that q ¢ 0 vanishes everywhere. We set V trap =0 for the region between the potential walls and we obtain the only remaining equation where μ is the chemical potential of the BEC and c 0 , the speed of sound in the BEC, is defined as Inserting the expansion y r d y into equation (4), and using equation (6) and all the properties of r 0 and θ 0 discussed above, we obtain the time-dependent BDG equations in first order in the perturbation δψ and its complex conjugate 5 Similar to the expansion of ψ, we can give an expansion of density and phase as r r a = + ( ) 1 0 and θ=θ 0 +f 0 +f. Via the equation y r = q e i , we can identify δψ and its complex conjugate with perturbations of density and phase as * a dy dy = + ( ) 2 and * f dy dy = -+ ( ) i 2 . From equations (8) and (9), we obtain In section 2, we found that the gravitational potential can be split into a stationary term and oscillating terms Here, we are only interested in the oscillating terms which lead to oscillating perturbations α and f. The stationary term leads to a stationary perturbation α c and f c , which slightly changes the evolution equations for α and f. This leads to effects of higher order in the gravitational potential. In appendix A, α c and f c are derived and evaluated for appropriate experimental parameters. In the following, they will be neglected and we will consider 2 . The first part of the last term in equations (10) and (11) vanishes inside the trap and becomes singular at the potential walls. Therefore, these terms can be neglected inside the trap and lead to the von Neumann boundary conditions a f ¢ = = ¢ 0 at the potential walls. We take the boundary conditions into account by using a mode decomposition 2 is a constant in first order in the perturbation, where  is the cross sectional area of the BEC in the y-z-plane. 5 Note, if we set Φ=0 and consider only points inside the trap, these equations lead to the time-independent BDG equations for i i that can be found in [24] and other standard books.
Projecting equations (10) and (11) on the mode j n (x) as As mentioned above, we only consider perturbations with wavelengths much larger than the healing length, which can be expressed as k n ζ=1. Therefore, the last term in equation (13) is much smaller that the second to last term and can be neglected. Then, by multiplying equation (13) by dV , we obtain that in first order in the potential perturbation. Taking this into account, taking the first time derivative of equation (13) and using equation (14) to replaceḟ n , we obtain for each n where we approximated because we assumed k l ζ = 1, and where we defined ω n =c 0 k n . This is a set of coupled, driven harmonic oscillators, which for vanishing driving, evolve with frequency ω n . Let us assume that the frequency of the external driving Ω is of the same order as ω n such that d w d Then, we obtain that the absolute value of the first term in the integral in equation (15)  . In this article, we will only consider situations in which just one mode l contributes to the evolution of the mode n via the coupling term in the integral in equation (15) 6 . Then, since k n ζ = 1 and k l ζ = 1, the second term in the integral in equation (15) will dominate significantly and we can neglect the first term. Finally, we find a set of ordinary, coupled linear differential equations for the amplitudes g n that are driven by the gravitational field there is direct driving through D n , the driving due to other excited modes through T nl and parametric driving through S n . Finally, we have to specify  dm f = -˙0. For that purpose, we project equation (10) on the constant function representing the density distribution of the ground state. It follows that  , which is the spatial average of the Newtonian potential. Then, we obtain for the driving moments, the parametric frequency modulation and the mode coupling coefficients respectively. We see the different signatures of the acceleration and the gravity gradient in the perturbation of the BEC. Via the direct driving D n , the gravity gradient couples to the symmetric modes while the acceleration couples to the anti-symmetric modes in the range of validity of our approximations. Similarly, acceleration couples modes of different parity and the gravity gradient couples modes of the same parity. Parametric driving S n affects all modes, but it is only due to the gravity gradient.
To obtain the full evolution of the phonon modes in the BEC, we have to include all significant damping effects. This is the content of the next section.

Damping
A phenomenological way to describe damping on the level of the GP equation was introduced in [26] and discussed further in [27]. We discuss this approach in appendix B. It is equivalent to adding a damping term gġ n n to the left-hand side of the harmonic oscillator equation (16), which leads to n n n n n n n l n nl l 2 In the following, we will discuss the two most important mechanisms of damping and their contribution to the damping constant γ n . There is a great deal of literature about damping of phonon modes in uniform BECs with periodic boundary conditions. In appendix D, we show that damping of phonon modes in a box potential with von Neumann boundary conditions can be described approximately using the expressions for damping in a uniform BEC with periodic boundary conditions. We will only discuss these expressions in the following.
For BECs at temperatures T such that k B T?ÿ ω n , where k B is the Boltzmann constant, one mechanism of damping is Landau damping, in which energy from the perturbations of the mean field is absorbed by the thermal bath of excitations. Landau damping in BECs was initially discussed in [28][29][30]. An expression for the damping constant γ n in a uniform BEC for general temperatures was derived in [31,32] and it is given as . For temperatures such that k B T?μ , we obtain the expression Note that (20) is linear in the temperature, while (21) is proportional to its fourth power, which means that the Landau damping can be lowered significantly by lowering the temperature further once the low temperature regime is reached. For the experimental parameters that we will consider in section 8 we have m  k T B and Landau damping is described by equation (21).
Another contribution to the damping rate that becomes most significant for low temperatures arises through Beliaev damping. The microscopic origin of Beliaev damping is the decay of a single phonon into two phonons of lower energy. The corresponding damping constant is given as [32]  is the Beliaev damping constant at zero temperature, and we assumed that the quotient of the atom density of the BEC and total atom density including the thermal cloud is close to one, which restricts our considerations to temperatures much smaller than the critical temperature of the BEC [22]. Note that (23) is proportional to the fifth power of the wave number while the Landau damping constant is linear in the wave number. This means that Beliaev damping becomes dominant for higher wave numbers. Since Beliaev damping remains even for zero temperature, we can conclude that the value of the wave number at which Beliaev damping becomes dominant is lower for lower temperatures.
The third mechanism of damping that has to be taken into account is damping due to atomic losses [35]; the atom-atom interactions and thermal fluctuations of the Bose gas lead to the evaporation of atoms from the condensate. The atom loss leads to damping and fluctuations of the phonon modes. In [35], it was shown that the numerical value of the corresponding damping rate γ loss is approximately the same as the evaporation rate of the condensate independently of the mode number. In general, the total damping rate can be written as

La Be loss
In the next section, we will discuss the number of phonons that correspond to the amplitude g n .

Average number of phonons in the mean field perturbations
The picture of perturbations of the order parameter ψ hides the true nature of the perturbations as phonons. However, the amplitude g n of the phase perturbation corresponds to a certain number of phononsN n,c that we can obtain by deriving the amplitude of the phase perturbation corresponding to a single phonon. For this purpose, we consider the energy of free fluctuations of the BEC density and phase without driving, which is given as where we used equation (5.76) of [23], the replacement J r a f = + ( ) i 0 and the equations of motion (10) and (11) for δV=0 and considering the von Neumann boundary conditions to be fulfilled. For the modes defined in equation (12), we obtain the energy , 2 6 n n a n n n is the number of atoms in the BEC and we neglected the contribution of the third term in equation (25), since ζk n =1. Without driving, equation (10) leads to . Furthermore, a solution of the equations of motion (16) without driving will be n n n for some phase j n . The same result is found for the steady state solutions on resonance later in equation (31). Therefore, using the linear dispersion ω n =c 0 k n , we find n n a n 2 2 2 The energy of a single phonon of mode n is given as ÿω n and we obtain the amplitude corresponding to a single phonon as z = ( ( )) g k N 2 2 n n a ,ph 1 2 . Additionally, we find for the average number of phonons in the coherent state created by the gravitational field of the moving mass pz In section 6, we will derive the amplitude of the dynamical modes of density perturbations induced by the moving mass and use equation (28) to extract the number of created phonons from the amplitudeḡ n . From equation (28), we can also find a condition on the maximal number of phonons in a mode for which our approach of considering δψ as a small perturbation is still suitable. We find that g n,0 10 −1 leads to  pz

Phonon creation for resonant driving
In sections 3 and 4, we derived the linearized differential equations for the evolution of the BEC under the gravitational influence of a moving source mass. In this section, we will derive the effect of the sinusoidal driving by solving equation (18) for particular cases.

Direct driving
The first case that we want to consider is the case of direct driving, when the initial excitation of all modes can be considered to be zero. Then, equation (18) reduces to g w j , 29 n n n n n n 2 where j j p = + 2 and The steady state solution on resonance for this driven and damped harmonic oscillator is if Ω is larger than ω n and j 0 =−π/2 for Ω smaller than ω n . We see that the steady state amplitude is F n ω n /γ n . From the solution (31), we see that there is a phase shift between the driving and the motion of the BEC of π/2 on resonance. Instead, if we consider times much shorter than the inverse damping rate, the solution is the undamped solution sin cos cos sin 32 n n n n n n , 0 2 2 n which, for resonance becomes cos sin cos . 33 n n n n n n 2 Hence, the amplitude grows linearly in time for . Together with equations (33) and (28) gives the number of coherent phonons created by the oscillating gravitational field of the massive sphere on resonance for times much shorter than the inverse damping rate of the mode under consideration as From the values for acceleration and gravity gradient in equations (3) and (34), we see that much less phonons are created due to the oscillating gravity gradient than through the oscillating acceleration when only the direct driving term D n is considered.

Mode coupling
Of course, in general, the solutions (33) and (31) are only reliable if the inter-mode coupling can be neglected. From the expression for T nl in (17), we see that the driving term due to a coherent excitation of mode l is of the same order as the direct driving term D n if If the inter-mode coupling cannot be neglected, either phonon pairs are created in modes n and l, or phonons in an excited mode will be shifted into a mode of higher energy. This can be seen by considering g l to oscillate with frequency ω l as l l l . Neglecting parametric driving and direct driving, we find the differential equation The right-hand side of equation (36) can be rewritten as Therefore, resonant driving of mode n is achieved if Ω=ω l +ω n or w w W = -| | n l , which correspond to the creation of phonon pairs and the shift of phonons between modes, respectively. These processes will appear as multi-mode squeezing and mode mixing, respectively, in section 7. Assuming initial excitation of mode n as n n ,0 0 , shorter times than the inverse damping constant and neglecting the back action on mode l, we find analogously to equation (33) on resonance for Ω=ω l −ω n . We see that, depending on the phase relation between driving and initial excitation of the mode l, the amplitude of the oscillation will increase from the start or it will decrease first and increase later. Hence, for an appropriately chosen phase relation, the mode coupling may be associated with a damping process. For values of  N 1 n lim,l , such damping processes may be stronger than the direct driving D n . In this situation, it may be more efficient when initial phonons are prepared in a certain mode and the gravitational field is measured through the induced loss of phonons from that mode. This possibility will be discussed in more detail in sections 7 and 8. We have to keep in mind that we neglected the back action on mode l, and therefore, equation (39) can only hold for short time scales.
For initially vanishing excitation, g n,0 =0, and j p = 2 l , we find for the number of created phonons in the mean field is the number of initial phonons in mode l that contribute to the mean field.

Parametric driving
As the third process of interest in this system, let us investigate the effect of the parametric driving term . For a classical system the parametric driving can only lead to an excitation of the system when the system is already excited. If we neglect all other processes and damping, start from an initial excitation that oscillates as for t<0 and consider parametric resonance Ω=2ω n , the amplitude after a time t is found as where g 0 is the initial amplitude. With equation (28), this leads to the expression for the average number of created phonons whereN n,0 is the initial number of phonons in mode n that contribute to the variation of the mean field.
In the next section, we will describe the phonons in the BEC as a quantum field and we will investigate the inter-mode coupling and the parametric driving in detail. We will recover equation (34) as the number of created coherent phonons, equation (40) from mode mixing and equation (41) from single-mode squeezing of a coherent initial state.

Quantum field description
For a description of the perturbations of the BEC as phonons, we start from a time-dependent external potential 7 s i n and . The external potential enters the Hamiltonian of the BEC as ò such that the full Hamiltonian is given by (see [23] for the free Hamiltonian) and  is the volume of the box potential. We neglect the contribution of the stationary part δV 0 as in section 3, consider the box trap potential V trap and the expansion 8 is the time-dependent energy shift of the ground state. Furthermore, we apply the Bogoliubov approximation, which replaces   a N a 0 , throughout the interaction process. Then, we obtain the grand canonical interaction Hamiltonian (see appendix Cfor the detailed derivation) The field operator is expanded as fulfill von Neumann boundary conditions at the potential walls (due to the vanishing density ρ 0 ) and are normalized with respect to the inner product and k n =πn/L were already defined before in section 3. This leads to ω n =c 0 k n for k n ζ=1.
sin , sin and the rotating wave 7 For other approaches to the description of BECs in time-dependent potentials see for example [36][37][38][39][40]. 8 Here we are using the absolute perturbation Ĵ in contrast to the description in section 3.
approximation, the interaction Hamiltonian can be written in the interaction picture as where the transition amplitudes are given as nn nn nn is reached for n=1. In the following, we will only consider processes on resonance; this means that n Ω ≔ LΩ/(πc 0 ) is an integer. Then, we find for the time evolution 2 , 2 a n d 2 . 5 8 n n n n n l n l n l n l n 0 , , Note that the second term in equation (57) only exists if n Ω is even. In the following, we will discuss all the different terms in equation (57), explain their meaning, investigate their effect on the phonon number and compare the results to those of section 3.

Coherent displacement
The first term in the exponent of the time evolution operator in equation (57) would create a coherent state when the other terms are neglected. This coherent state leads to non-vanishing expectation values of the quadratures 2 , which can be identified with r a 0 and r f 0 , respectively. Therefore, the first term in equation (57) can be identified with the direct driving D n . Taking only this term into account, the average number of phonons is given as  . Hence, for a highly excited initial state, a small amplitude M nn still can lead to a measurable change in the average number of phonons. Such initial phonons may be created by applying an external electromagnetic linear or harmonic potential to the BEC superimposed with the trap potential by the same mechanisms that we consider here for the measurement process.

Multi-mode squeezing
Additionally to one mode squeezing, we find multi-mode squeezing with the third term in the exponent of the time evolution operator in equation (57). In particular, this process can be induced by acceleration and the gravity gradient, while acceleration only creates phonon pairs in modes with different parity, the gravity gradient creates phonons in modes with the same parity. For small mode numbers n and l and n+l even, the amplitude M nl will be of the same order as M nn . However, the quotient of the amplitude M nl /M nn increases with increasing mode numbers l and n when n−l is kept constant. Therefore, to generate phonons with higher frequencies, multi-mode squeezing will be more efficient than one mode squeezing. In all cases, the resonance condition Ω≈ω n +ω l has to be fulfilled to create phonons. Only taking multi-mode squeezing into account, the evolution of the creation and annihilation operators can be represented as using the twomode squeeze operator T l n n l n l l n l l n l 2 , Let us assume that the system starts in an initial displaced squeezed state a x ñ | , n n,0 with n<n Ω . Only taking multi-mode squeezing into account, we obtain the state a x ñ ( )| S t , T n n,0 . We have = +- T n T nn n n n n n n n n , , T n n T nn n n n n n n n n , , and the total average number of phonons in mode n and n Ω −n becomes

Mode mixing
The last term in the exponent of the time evolution operator in equation (57) leads to mode mixing, which partially corresponds to the mode coupling process; the second term on the right-hand side of equation (16).
The amplitude is the same as for the two-mode squeezing and the process will be efficient for high mode numbers n and l when n−l is small. When only the mode mixing is taken into account and only resonant processes are considered, the time evolution can be described by a beam splitting operation l n l n l l n l l n l n l

,
This operator couples all modes with a distance n Ω . This means, that we obtain n Ω −1 systems of coupled modes 9 . Starting with a displaced squeezed initial state a x ñ | , n n,0 for a single-mode n and assuming that ,c,0 is the total number of phonons in the initial state.

Quantum Cramér-Rao bound
From our analysis above we see that higher creation rates of phonons are obtained for initially excited states. From the created phonons the gravitational field amplitudes a Ω and W G can be inferred. Therefore, optimal precision for the measurement of the gravitational field via the creation of phonons is obtained for excited states, which we can call probe states. The estimation via the squeezing effect represented by M nn and M ln is sometimes called a squeezing channel. The estimation via the mode mixing represented by A ln and B ln is called a mode mixing channel. It was shown in [45] that the optimal estimation of a parameter via either of the two channels is obtained for a probe state where all particles are squeezed particles. Here, optimal means highest precision per phonon in the probe state. In the following, we will give this optimal precision in terms of the quantum Cramér-Rao bound, which gives an upper limit for the precision of the estimation of a parameter via an estimation channel for a specific initial state for all possible measurements. We obtain the limit for the absolute precision for a measurement of the driving parameter  where # rep is the number of consecutive independent measurements and I ò is the quantum Fisher information for the driving parameter ò.
Let us assume that we start with two non-entangled modes each in a squeezed probe state with squeezing parameters r 0,1 and r 0,2 . The quantum Fisher information for the optimal estimation of the driving parameter ò via the two-mode squeezing channel is given as   [45]. Both expressions correspond to the Heisenberg limit and represent an upper bound for the achievable sensitivity if the probe state is set up as described above.
In the next section, we will investigate the sensitivity for the estimation of acceleration and gravity gradient for examples of experimental parameters.

Experimental parameters and measurement sensitivity
In this section, we give necessary experimental parameters for the measurement of the gravitational field of an oscillating mass via phonon creation in a BEC. To provide an example, let us consider rubidium BECs and ytterbium BECs. For Rb-87, we have an interaction constant of l »´-1.3 10 m Rb 7 (calculated from the measured scattering length a scatt =λ/8π≈98a 0 reported in [46], where a 0 is the Bohr radius). The mass of a Rb-87 atom is 1.44×10 −25 kg. The scattering length of Yb-168 can be found in [47] and leads to a interaction constant of approximately λ Yb =3.35×10 −7 . The mass of a Yb-168 atom is 2.79×10 −25 kg. In the following, we will assume that the cross section of the BEC in the y-z-plane is circular.
The experimental time t exp is limited by the half-life of the BEC and the half-life of the phonons. We assume that the half-life of the phonons, which is proportional to the inverse damping rate, is much larger than the halflife of the BEC. We find that this assumption is met for the experimental parameters chosen. Therefore, the halflife is the only limit for t exp in the following consideration. In section 5.4 of [24] and in [48,49] it is shown that where D is the decay constant. This implies a quadratic dependence of the half-life on the inverse density. More precisely r = ( ) t D 3 2 hl 2 . For example, in an experiment with rubidium atoms [50], the corresponding decay constant was found to be 1.8×10 −29 cm 6 s −1 . In an experiment [51] with ytterbium an even smaller decay constant of 4×10 −30 cm 6 s −1 was found. Therefore, we assume that the density is less or 9 These systems of coupled modes can be seen as quantum Markov chains that are infinite on one side. B(t) gives rise to the time evolution of an initial state on this Markov chain. The beam splitter operation is a completely positive map [43] and Markov chains can be understood as completely positive maps between sites of a graph [44]. Here we have a one-dimensional undirected graph with one end. equal to 10 13 cm −3 in order to achieve a half-life of at least t hl =100 s. This allows the assumption that the number of atoms can be considered to be constant during a single run of the experiment for an interaction time of t exp =10 s. As we discussed in section 4, the decay of the BEC leads to a damping of the phonons inside the BEC, with a value for the damping rate g loss that approximately matches the value of the BEC decay rate. From the discussion above, we find that γ loss ≈10 −2 s −1 is a conservative assumption. For all the experimental parameters that we will consider in the following, we find that the Landau damping rate g n L is of the order of --10 s 3 1 and Beliaev damping is significantly smaller. Therefore, the total phonon damping rate is dominated by the atom loss, we can set γ n =γ loss , and we obtain g - t n 1 exp . We assume that  p p W´=t 2 5 2 0.5Hz exp , which means that at least about 5 cycles of the driving oscillation elapse during one run of the experiment. To fulfill the condition k B T = μ, we consider a temperature of 1 nK. In a uniform BEC this leads to a relative depletion 10 of the density of atoms in the ground state of the order of 10 −4 . Accordingly, less than -N 10 a 3 atoms are in thermally excited states. The back action of these atoms on the condensate leads to a stationary deformation of the order parameter and the number of phonons created due to the oscillating gravitational field is changed slightly. Hence, thermal depletion is a small effect for the parameters considered in this article that can be neglected.

Phonon creation due to direct driving
In the following, we want to find experimental parameters that allow for the detection of gravitational acceleration by measurement of phonons created via direct driving. To achieve the maximal creation rate, the mode number n has to be as small as possible. Therefore, we consider n=1 for the creation of phonons by direct driving due to the oscillating acceleration. From the expression for the frequency ω n =c 0 k n and the definition of the speed of sound in equation (7), we find The resonance condition is ω 1 =Ω and we assumed that Ω2π×0.5 Hz. This condition is fulfilled if we choose L=200 μm, which corresponds to ω 1 =2π×1.5 Hz for rubidium and ω 1 =2π×1.2 Hz for ytterbium. For a successful detection, we need the signal to noise ratio  = D  [52], the temperature of the BEC varies by about 20% between two runs. We assume that the temperature is about 1 nK, which implies that k B T?ÿω n . Hence, the Bose-Einstein statistics tells us that the value for the variance of the number of thermal phonons is approximately equivalent to the value for the average number of phonons. We . We will discuss possibilities of a detection process in section 9. Furthermore, we assume that the number of atoms varies by about 10% between two experiments, which means DÑ N 0.1 a a . Since each experiment takes about 10 s, we can consider a number of repetitions of # = 10 rep 4 , which corresponds to two days of consecutive measurements. For example for ω 1 ∼2π×1 Hz, we obtain thatÑ 20 n,th . Then, the thermal fluctuations are the main source of fluctuations, and we find that the creation of a single phonon by the gravitational field would be sufficient to reach a signal to noise ratio of the order of 10. In the following, we assume  = 10 SNR,1 , and we give the experimental parameters necessary to achieve this goal. The necessary number of created phonons for the detection processN n cr , can be used to give a lower bound for the number of atoms in the BEC by using equations (70) and (34), which leads to 10 For the depletion, we used the formula (4.50) in [23]: , where ρ 0 (T) and ρ 0 (T=0) are the density of atoms in the ground state at temperature T and T=0, respectively. The fixed length and the fixed number of atoms can be used to fix the ratio between the transversal diameter d and the length of the BEC. We obtain The density should be as large as possible in order to keep d small. Therefore, we match the upper bound for the density given above and set ρ=10 13 cm −3 .
In figure 2, a table can be found in which some example values for experimental parameters are listed that could be used for the detection of the gravitational field of an oscillating gold or tungsten sphere of mass M for the two cases of M=200 g and M=0.2 g. Values for the parameters are given for both a rubidium and an ytterbium BEC. For a source mass of 200 g, we find that the phonon creation should be observable with state of the art technology. In the case of ytterbium, the number of phonons for which the mode coupling becomes significant is smaller than the number of created phonons. This suggests that a regime can be reached in which mode coupling and parametric driving may supply an alternative detection scheme. We will investigate this possibility in the next subsection. If we repeat the above calculations for the creation of phonons due to the oscillating gravity gradient, we find that we would need about 10 10 atoms for the detection using the direct driving process. This is experimentally out of reach at the moment. Mode coupling and parametric driving will be much more efficient for gradiometry.

Phonon creation due to squeezing
In section 7, we discussed the creation of coherent phonons in an initial coherent state and the creation of additional squeezed phonons in an initial squeezed state due to the squeezing processes. The first situation corresponds to the creation of phonons in the mean field due to the parametric driving of the BEC by the oscillating gravitational field discussed in section 6. We found that the total average number of phonons that are created on resonance in the mode n by multi-mode squeezing is given as n n n n n n , , , whereN n,0 was the initial number of phonons (squeezed plus coherent) in mode n and n Ω =Ω L/(π c 0 ) is an integer since we consider resonant driving.
The amplitude and ζ as small as possible. In contrast to the direct driving, the number of atoms in the BEC does not appear in the number of created phonons due to squeezing. Since z l r = 1 0 , the healing length can be minimized by choosing for a density of ρ 0 =10 13 cm −3 , which we identified as the maximum when an experimental time of at least 10 s is to be obtained. Note that only gravitational acceleration contributes for n Ω odd and only the gravity gradient contributes for n Ω even. We can minimize --W | ( )| n n n 2 2 by choosing n=(n Ω +1)/2 if n Ω is odd and n=n Ω /2+1 if n Ω is even. We consider n=2 and n Ω =3 for the measurement of the acceleration and n=3 and n Ω =4 for the measurement of the gravity gradient in the following. For the detection of the oscillating acceleration, we fix the length of the BEC to m = L 200 m. For the detection of the oscillating gravity gradient, we have to increase the length to obtain more realistic values for the other experimental parameters; we consider m = L 500 m. Figure 2. COHERENT/DIRECT DRIVING: this table shows some generic values for the experimental parameters necessary to detect phonons in a BEC created by direct driving due to the oscillating gravitational acceleration of amplitude a Ω induced by a small oscillating sphere of mass M with a signal to noise ratio of the order 10. It is assumed that the density of the BEC is r = -10 cm 13 3 , the length of the uniform trap potential well is m = L 200 m, the temperature of the BEC is T=1 nK and the measurement precision is of the order of a single phonon. The interaction time for each experiment is assumed to be t exp =10 s. In each case, about 10 4 repetitions of the experiment are considered. Under these conditions, the minimal distance R min to the BEC, the oscillation amplitude δ R , the frequency of the driving Ω/2π, the number of atoms N a in the BEC, the ratio of length of the BEC and healing length L/ζ and the ratio of the length and the diameter d/L of the trap potential are given. Additionally, the table shows the average number of created phonons, the value for the number of phonons at which the inter-mode coupling becomes significant N 2 lim,1 and the average number of thermal phonons. Now, we can calculate the number of initial phononsN n,0 that are necessary to createN n cr , phonons on average in a single experiment, which we need to achieve a signal to noise ratio of the order of 10 after 10 4 repetitions of the experiment. We The relation between the signal to noise ratio andN n cr , can be derived from equation (40) as The results for the above parameters are presented in the table in figure 3. Additionally to multi-mode squeezing, we can also consider the utility of single-mode squeezing for the detection of the gravitational field of the oscillating mass. From equation (54), we see that single-mode squeezing can only be induced by an oscillating gravity gradient. The resulting parameters for single-mode squeezing can be found in figure 4. We see that single-mode squeezing gives better parameters than multi-mode squeezing for the measurement of the gravity gradient. This is because, for small W ( )

Measurement via mode mixing
The remaining channel that can be used for a measurement is mode mixing; phonons in one mode will be transferred to another mode due to the oscillating gravitational field. Since the mode frequencies are equidistant, if there is a driving frequency w w W = - . Therefore, phonons that are transferred to a higher mode from an initially excited mode will not stay in that mode but will be transferred up the whole cascade of resonant modes. Hence, starting from an excited state in Figure 3. TWO-MODE SQUEEZING: these tables show some generic values for the experimental parameters necessary to detect the phonons in a BEC created by two-mode squeezing due to the oscillating gravitational acceleration of amplitude a Ω and the oscillating gravity gradient of amplitude W G induced by small oscillating sphere of mass M with a signal to noise ratio of the order of 10. It is assumed that the density of the BEC is r = -10 cm 13 3 , the temperature of the BEC is T=1 nK and the measurement precision is of the order of a single phonon. Furthermore, the length of the uniform trap potential well is L=200 μm for the consideration of a Ω and m = L 500 m for the consideration of W G . The interaction time for each experiment is assumed to be t exp =10 s. In each case, about 10 4 repetitions of the experiment are considered. Under these conditions, the minimal distance R min to the BEC, the oscillation amplitude δ R , the frequency of the driving Ω/2π and the minimal number of atoms N a min in the BEC are given. Additionally, the tables show the values for the squeezing parameter r 1,2 and r 1,3 , the average number of created phonons N 2,cr and N 3,cr , the necessary number of initial phononsN 2,0 andN 3,0 and the average number of thermal phonons. one mode, we could measure the decrease of the number of phonons in this mode or the increase of the total number of phonons in all modes.
From the amplitude in equation (54), we see that we can apply the same arguments as for the multi-mode squeezing, the only difference being that we consider the two modes n and -W n n . If we consider n=3 and n Ω =1 for measurement of the gravitational acceleration, we obtain less favorable experimental parameters than those in the table in figure 3 above. If we consider n=3, n Ω =2 for the measurement of the gravity gradient, we recover the same parameters as in the table in figure 3 with the exception of the driving frequency which is decreased.

Quantum Cramér-Rao bound
We can obtain an upper bound for the sensitivity of the measurement of oscillating gravitational fields using phonons in BECs by considering the quantum Cramér-Rao bound that we introduced in equation (69).

Conclusions and discussions
The necessary experimental parameters for the measurement of the gravitational field of an oscillating sphere of mass M=200 g due to the direct driving of phonon modes in a BEC with a signal to noise ratio of the order of 10 seem to be ambitious but not inaccessible (see the first table in figure 2 for details). State of the art experiments with ultracold rubidium BECs (at about 1 nK) use a number of atoms of the order of 10 5 [1,21,53,54] and atom numbers of the order of 10 6 are planned for a new generation of experiments [55][56][57]. In section 8, we argued that the interaction time for a single experiment of the order of 10s can be achieved by choosing a low atom density of the order of 10 13 cm −3 . The parameters for the case of M=0.2 mg are out of range of state of the art experiments; the number of 10 8 atoms necessary to achieve detection, with a signal to noise ratio of 10, is not obtainable. Nevertheless, this parameters may be achievable in the future.
Besides phonon creation due to direct driving, we investigated phonon creation due to parametric driving resulting in squeezing and mode mixing. This driving mechanism turns out to be of advantage when the phonon modes are initially in an excited state. Such initial excitations may be created by adding an oscillating external electromagnetic linear or harmonic potential to the already existing BEC trap potential. Then, phonons would be created by the mechanisms that we consider here for the measurement process. For example, this could be the direct driving or parametric driving from the vacuum, where the latter is equivalent to the dynamical Casimir effect in Bose-Einstein condensates [42]. See also [58], where parametric amplification of excitations of phonons modes due to a modulation of the transverse trapping frequency of a BEC is discussed in detail.
We gave necessary experimental parameters for the measurement of the gravitational field of an oscillating massive sphere using parametric driving in figure 3. Firstly, it is interesting to note that our theoretical considerations predict that ytterbium would perform much better than rubidium. Secondly, even for a few initial phonons, the parametric driving mechanism is more efficient than the direct driving as a lot less atoms are needed in the BEC to achieve a signal to noise ratio of the order of 10. This is particularly useful for the measurement of the gravitational field of smaller masses; the acceleration due to a sphere of 200 mg can be G induced by a small oscillating sphere of mass M with a signal to noise ratio of the order of 10. It is assumed that the density of the BEC is ρ=10 13 cm −3 , the temperature of the BEC is T=1 nK and the measurement precision is of the order of a single phonon. Furthermore, the length of the uniform trap potential well is m = L 500 m. The interaction time for the experiments is assumed to be t exp =10 s and about 10 4 repetitions of the experiment are considered. Under these conditions, the minimal distance R min to the BEC, the oscillation amplitude δ R , the frequency of the driving Ω/2π and the minimal number of atoms N a min in the BEC are given. Additionally, the tables show the values for the squeezing parameter r 1 , the average number of created phonons N 1,cr , the necessary number of initial phononsN 1,0 and the average number of thermal phonons. measured with 10 6 atoms when the initial state contains about 100 coherent phonons. The measurement of the gravity gradient of small spheres with masses of the order of 200 g or less using direct driving is completely out of reach. However, it can be achieved with parametric driving with a BEC consisting of about 10 6 atoms and 10-100 initial phonons. In this article, we assumed that the BEC is always much smaller than the distance between its center and the center of the source mass. It would be interesting to relax this condition in a future investigation. On the one hand, decreasing the distance between the source and the BEC beyond that limit may lead to further improvement in the measurement sensitivity. On the other hand, measurements close to the surface of small masses may allow for experiments to search for hypothetical fifth forces or to measure Casimir-Polder forces. Finally, we calculated the quantum Cramér-Rao bound for the measurement precision considering the same experimental parameters. Hence, there seems to be a lot of potential for improvement by taking the insights of quantum metrology into account. This potential may become accessible through measurement schemes other than just the direct counting of created phonons such as, for example, a homodyne measurement of phonon modes.
In the this article, we assumed that the measurement technique employed reaches single phonon sensitivity. We are optimistic that single phonon sensitivity can be achieved in the near future if research effort is made in this direction. In particular, it seems that a precision of tens of phonons is achieved in experiments like the ones presented in [59,60]. Single phonon sensitivity may be achievable as experimental procedures evolve. There is a variety of possibilities to measure phonons in BECs, some of them may, in the future, give a precision high enough for our purposes. In particular, one can either try to measure the phase or one can try to measure the density, which are conjugate variables and contain the same information. An interesting approach for measuring the phase is presented in [61]. It is denoted as 'heterodyne detection' by the authors: after the modes are excited, the trap potential is switched off and the BEC starts to expand and fall freely. During the expansion the energy contained in the phonons is transformed into the kinetic energy of atoms. These free particles interfere with the atoms in the ground state. The interference fringes contain the information about the phonons. Numerical simulations and an approximate analytical derivation for this process are given in [62].
Another option for a detection scheme would be time of flight measurements, where phonons are mapped to horizontal atomic momenta and, after a certain time of vertical free fall in the gravitational field of the Earth, the momenta can be read from the horizontal position of the atoms. A third option for the measurement of phonons in BECs would be direct light phonon couplings. For example in [63], stroboscopic measurements of phonon modes were considered for the creation of squeezing and entanglement of phonon modes. In [64,65] nondestructive phase-contrast imaging was used to observe the bulk perturbations of a BEC. Finally, a fourth option for the measurement of phonons would be the coupling to atomic quantum dots submersed in the BEC [66]. It would be very interesting if experiments could be performed to investigate the sensitivity of different measurement schemes for phonons. A first step towards an experimental realization of our proposal could be experiments using one of the above techniques to simply measure the thermal spectrum of phonons in a BEC with high precision. A second step could be to create phonons in the BEC by Bragg scattering of laser pulses or by periodic modulations of the trapping potential and try to measure them on top of the thermal spectrum. In a last step, the interaction with an oscillating source mass can be implemented.
Additional noise sources that we did not discuss in the main part of this article are Newtonian and seismic noise that give rise to acceleration noise a x noise . The Newtonian noise also introduces a noise term into the gravity gradient, which we assume to be negligible since the sources of gravity gradient noise will be far away in comparison to the extension of the BEC, which means that the gravity gradient noise will be highly suppressed in comparison to the gravitational acceleration noise 11 . A generic example of the square root of the displacement spectral density, S x 1 2 , in a modern laboratory environment close to traffic is shown in figure 3.3 of Tobias Westphal's PhD Thesis [67] for the case of the physics department of the university of Hannover. The square root of the displacement spectral density S x 1 2 for 1 Hz is of the order --10 m Hz 7 12 12 . We can assume that the laboratory structure is not driven resonantly in this frequency range and that damping can be neglected. Then, the susceptibility can be approximated as 1/ω 2 , and we find an acceleration spectral density of about w =~- . Hence, the Newtonian and seismic noise background has to be lower by only about one order of magnitude to get below the order of the gravitational acceleration due to a 200 g source mass. For the case of M=200 mg, the Newtonian and seismic noise have to suppressed by two orders. Both situations should be achievable by choosing a quieter environment, e.g. a site underground far from human induced noise, and a vibration isolation chain [68]. As an example, advanced LIGO is engineered to achieve noise levels at the mirrors of about´--5 10 m Hz 19 1 2 [69] for frequencies above 10 Hz, which would 11 This is a clear advantage of measurements of the gravity gradient of small objects besides the property that the gravity gradient of a sphere close to its surface is independent of its radius and only depends on its mass. 12 Actually, there is a dip slightly below 1 Hz, which the experimental parameters may be tuned into to lower the seismic noise background. be more than sufficient to make the acceleration noise negligible in comparison to the gravitational acceleration induced by the source mass. In particular, vibration isolation chains have to be included in any case since the source and the detector must not be coupled significantly through the devices holding them at their respective positions.  . We assume that the number of atoms in the BEC is of the order 4×10 6 , the length of the BEC is m = L 300 m and its density is r » . We find that f q~-˙10 c 0 8 . Therefore, all terms proportional to α c and f c that can appear in equations (10) and (11) are negligible, and we are justified to treat the effect of the sinusoidally time-dependent terms in Φ independently of α c and f c .
The basic idea of the phenomenological treatment of dissipation presented in [26,27]  0 for a given chemical potential μ. Therefore, the time evolution operator for the damped system out of equilibrium has to vanish identically on this state. This operator is derived by removing the chemical potential from the differential operator  y [| | ] 0 2 and multiplying the resulting differential operator with the factor (1+iΛ), where Λ is the damping constant. Then, the resulting non-unitary differential operator vanishes on the equilibrium state as wanted.
From the GP equation, we obtain the time-independent GP equation when iÿ∂ t is replaced by μ. The equilibrium state is the initial state y r = Following the steps in section 3, we obtain equation (18) in first order in Λ.
ij n n i j i j i j n i j i j i j 3 0 * * * * * * * * * * * * ò f i i i , , x y z and so on. We are considering interaction times that are only a factor 5 larger than the inverse frequency of the modes under consideration in this article. Therefore, the delta function in equation (D1) have to be replaced by a sinc-function as [25] The width of the sinc-function w w w + + ( ( ) ) t sinc 2 i j exp 0 is still small enough to justify that we can neglect the third term in equation (D1). Let us assume that the BEC trap has a square cross section of edge length L tr . In analogy to the expressions for the phonon mode functions for a BEC in a box in the x-direction given in equation (50) The sinc-functions in equation (D1) only deliver significant contributions in combination with M n ij when the momentum relation agrees with the energy relation in the argument of the sinc-function. In  For n z of the order of 1, the relation i z +j z =n z in combination with i x =j x or i y =j y will lead to a value for ω n −ω i −ω j much larger than the width of the sinc-function. Therefore, we can write We find that the Landau damping rate for phonons in a BEC in a box potential is of the same order as the Landau damping rate in a uniform BEC with periodic boundary conditions. For the Beliaev damping, we find an increase of about one order. Since Beliaev damping is strongly suppressed for our parameter range, we conclude that we can use the expressions for the damping constants derived for periodic boundary conditions to describe the BEC in a box potential.