Analogue simulation of gravitational waves in a 3+1 dimensional Bose-Einstein condensate

The recent detections of gravitational waves (GWs) by the LIGO and Virgo collaborations have opened the field of GW astronomy, intensifying interest in GWs and other possible detectors sensitive in different frequency ranges. Although strong GW producing events are rare and currently unpredictable, GWs can in principle be simulated in analogue systems at will in the lab. Simulation of GWs in a manifestly quantum system would allow for the study of the interaction of quantum phenomena with GWs. Such predicted interaction is exploited in a recently proposed Bose-Einstein condensate (BEC) based GW detector. In this paper, we show how to manipulate a BEC to mimic the effect of a passing GW. By simultaneously varying the external potential applied to the BEC, and an external magnetic field near a Feshbach resonance, we show that the resulting change in speed of sound can directly reproduce a GW metric. We also show how to simulate a metric used in the recently proposed BEC based GW detector, to provide an environment for testing the proposed metrology scheme of the detector. Explicit expressions for simulations of various GW sources are given. This result is also useful to generally test the interaction of quantum phenomena with GWs in a curved spacetime analogue experiment.


I. INTRODUCTION
In the 36 years since the seminal proposal of Unruh to measure an acoustic analogue to Hawking radiation from a "sonic horizon" in a fluid [1], interest in analogue simulation of gravitational fields has grown from theoretical proposals to experiments in numerous systems. These include Bose-Einstein condensates (BECs) [2][3][4], water waves [5,6] and optical fibres [7] among others. Particular interest has been shown in using the phonon field in a BEC, since this is a quantum system and so allows for the study of how non-classical properties, such as entanglement, are modified or generated by simulated gravitational fields.
This has recently culminated in the first observation of the entanglement of acoustic Hawking radiation [8], potentially providing clues to fundamental questions for quantum gravity, such as the information paradox. In addition to Hawking radiation from a waterfall horizon, other proposed simulations using BECs have included conformal Schwarzschild black holes [9,10], rotating black holes [11], FRW geometries [12,13], inflation [14,15] and extensions to Einstein's general relativity, such as aether fields [9]. In past work [16], two of us have considered the simulation of gravitational waves (GWs) in 1 + 1 dimensions. These are perturbations of spacetime generated by a changing quadrupole moment of a mass distribution, and have recently been detected in a milestone moment in science [17][18][19][20][21]. This has led to a new field of GW astronomy, enabling the exploration of the Universe through gravitational as well as electromagnetic radiation. Simulating GWs in fluid systems could be of astronomical interest, for example, in studying GWs in numerically challenging, strong-field regimes.
Furthermore, since BECs are quantum systems, this could enable the study of predicted effects such as particle creation in GW backgrounds [23], and quantum decoherence due to GWs [24]. While we are interested in simulating the effect of GWs, simulating the evolution of a GW itself on a curved background has also been proposed in [22], where a metamaterial emulates a curved background space-time and two-photon states model the evolution of a GW on that background.
Simulating GWs with BECs could also be used in studies of a proposed BEC GW detector [25][26][27]. This detector consists of a BEC constrained to a rigid trap with a prepared quantum state of phonons, such as a two-mode squeezed state. The transformation induced by the GW produces mode-mixing and phonon creation in a phenomenon resembling the dynamical Casimir effect [28,29], making the final state distinguishable from the probe state, i.e. decreasing the fidelity between the initial and final states. The lower the fidelity between the probe and final state, the better the estimation. Non-classical squeezed states allow for quantum metrology techniques resulting in better estimation than a classical device. At resonance, mode-mixing or phonon creation are maximised giving rise to optimal parameter estimation. Such a quantum resonance process is absent in laser interferometers since the frequencies of the GWs are far from the optical regime. However, using resonance to detect GWs was the concept behind the first GW detector proposals, Weber bars, which are metal objects measuring metres in length. The GW resonance in the BEC detector is similar since the much smaller size, O (µm), is compensated for by a much smaller speed of sound, O (mm/s) compared to O (km/s) [30]. However, the BEC detector can be cooled to considerably lower temperatures, O (nK), and is a strictly quantum device. This allows for the use of quantum metrology and, therefore, sensitivities that are inaccessible to classical devices [31]. A discussion of the viability of such a detector and further details can be found in [25]. Further studies of the viability of such a detector are in progress and this article is part of this effort. Simulating the effect of a GW derived in [25] could be useful for testing the metrological scheme proposed in [25,26].
Here we extend the work on 1 + 1 GWs to the simulation of 3 + 1 GWs in BECs in a covariant formalism [16]. With this extension, all properties of a GW can be simulated, such as its polarization and propagation vector, and the conformal factor in front of the analogue metric no longer diverges. The paper is outlined as follows: in Section II we present the spacetime metric of a GW and the acoustic metric in a BEC. In Section III we derive and demonstrate the simulation of a GW metric in 3 + 1 dimensions, as well as the metric derived in [25]. Examples of GWs are presented in Section IV, giving explicit forms of flow velocities needed to simulate the effect of commonly investigated GW sources, including compact binary inspirals and neutron star spin down. Section V reduces the metric derived in Section III to 1 + 1 dimensions, and compares this to previously published work in [16], and we conclude in Section VI.

A. Definitions and conventions
Throughout this paper, we use the metric signature (−, +, +, +), the coordinates used are Minkowski coordinates given by (ct, x, y, z) unless otherwise stated, and the Minkowski metric in these coordinates is given by (1)

II. GRAVITATIONAL WAVES AND THE ACOUSTIC METRIC
A. GW spacetime metric We first consider the general form of the metric tensor perturbed by a GW. For a single source GW far from the source, the metric tensor can be expressed as [32,33] g (gw) where η µν is the flat Minkowski metric defined in Section I A, and h µν is some perturbation corresponding to the passing GW, parameterised by ǫ, where |ǫ| ≪ 1. Standard notation omits this ǫ and applies the condition |h µν | ≪ 1, but we use ǫ here as a global perturbation scale factor for consistency and clarity.

Transverse traceless gauge
This perturbation h µν can be expressed in the transverse traceless (TT) gauge, in coordinates x µ T T , for a GW travelling in theẑ direction as [33] where h + and h × are time-dependent functions corresponding to two "polarisations" of the GW. These h + and h × functions are typically called "strain" functions. We ignore the z dependence of the strain functions, as the wavelength of a GW is typically much longer than the width of a BEC (for example, the GWs detected by the LIGO and Virgo collaborations have wavelengths exceeding 10 6 m). Outside the source of the GWs, these strain functions obey the simple wave equation We introduce the GW metric in this gauge, as it is the clearest and most widely known, despite not necessarily being the most physically useful.

Fermi normal coordinates
Any metric in linearised gravity of the form of Eq. (2) has a "gauge freedom", namely a choice of small coordinate transformation in any arbitrary direction. Consider a linearised coordinate transformation with some function ζ, such that Under such a coordinate transformation, the metric in Eq. (2) transforms as Hence, making a coordinate transformation x µ T T → x µ = x µ T T + ǫζ µ from the TT gauge with the function the metric perturbation in Eq. (3) in these new coordinates is where I n is the n dimensional identity matrix, and These coordinates are Fermi normal coordinates, and are the inertial frame limit of the "proper detector frame". This metric perturbation can also be derived by considering coordinates matching proper length and time, then linearising the metric with respect to

B. Acoustic metric
To simulate a GW metric in a BEC, we will follow the description of a BEC on a general background metric given in [13,34,35]. This description models the BEC as a barotropic, irrotational and inviscid fluid, in a covariant formalism. We are interested in i) simulating a spacetime metric using a quantum system and ii) simulating the effects of spacetime dynamics on a phononic field. In both cases we require a covariant formalism that enables us to properly describe a general relativistic spacetime and its effects on quantum fields. The formalism developed in [13,34,35] enables us to do so. We point out that the system that we consider here is a regular BEC, as those currently demonstrated in the laboratory. This system is usually described with non-relativistic quantum mechanics. However, in ii) we are taking into account the underlying spacetime background on which the BEC sits on, which requires the covariant treatment mentioned above. We point out that we are not considering a system that is moving with relativistic speeds or has excitations with relativistic energies. Such a relativistic system would also need to be described by the same formalism since covariance is also necessary.
Any BEC which can be described as a superfluid is automatically barotropic and inviscid.
In the superfluid regime, the BEC is described by a classical mean field φ expressed as with quantum fluctuationsψ, defined in terms of the total fieldΦ aŝ Φ = φ 1 +ψ .
We are interested in the behaviour of these fluctuations in the "phononic" regime. The relativistic phononic regime condition can be written explicitly as [13] |k| where k is the spatial frequency of a phononic excitation of the BEC, u is the flow velocity of the BEC defined as and ξ is the "healing length" defined as λ encodes the strength of the interaction, defined in terms of the interaction potential U as where extra terms are φ 6 interactions and higher, which we ignore here. The interaction strength λ is related to the scattering length a by Taking the non-relativistic limit of this condition, we find that the phonons should have wavelengths far longer than the healing length ξ. In this regime, and with certain additional assumptions about the mean field properties, the fluctuations obey a relativistic Klein- for a tensor G µν with determinant G, called the "acoustic metric". The general acoustic metric is given by (see Appendix B) where g µν is the background spacetime metric, ρ is the bulk density defined as ρ = φ * φ, r is related to the speed of sound c s as and v is the normalised flow velocity defined as The speed of sound c s is defined as where The flow velocity v is normalised as Note that the definition of the flow velocity imposes irrotationality, i.e.
Due to the global phase symmetry of the Lagrangian that Eq. (18) is derived from, there is a conserved current. The conservation of this current can be expressed as also called the continuity equation, where ∇ µ is the covariant derivative with respect to g µν .
The velocity normalisation |u| and density ρ can also be directly related to the internal and external potentials U and V as

III. GRAVITATIONAL WAVE SIMULATION
In this paper, we present two results corresponding to different types of simulation. The first in Section III A is a direct simulation of the GW metric in Eq. (2). The second result in Section III B is a simulation of the acoustic metric derived in [25], to test the proposed metrological scheme.

A. GW metric simulation
The goal of this section is to directly simulate a GW spacetime, such that the acoustic metric G µν has the form We start with the GW metric in Fermi normal coordinates, as introduced in Section II A 2.
If we consider a BEC at rest in these coordinates, i.e. v µ = −cδ 0 µ , where the background metric is the flat Minkowski metric, then the acoustic metric has the form It should be noted that the density is not completely unrestricted; the choice of flow velocity restricts the density through the continuity equation. From the definition of the flow velocity in Eq. (14) and the choice of normalised velocity above, which implies that |u| can only be a function of time. Hence, with the continuity equation (Eq. (B16)), for this particular choice of normalised velocity, we must have

Bulk properties
Comparison of Eq. (8) and Eq. (29) suggests that, to simulate the GW metric, we must modulate the speed of sound as where c s0 is the speed of sound in the absence of a simulated GW. If we rescale the time coordinate as then the acoustic metric has the form which matches the desired metric in Eq. (28) up to a conformal factor. The conformal factor will be discussed further later in this subsection, as well as in Section III C.

Implementation
As mentioned earlier in Section III A, we must have |u| dependent on time only. Furthermore, having a density that is changing in time in the absence of flows implies changing the local atom number density in the BEC in a uniform and precisely controlled way, which seems experimentally unfeasible. If it is possible to control the density independently of the flows, then modulating the density to match the speed of sound sets the conformal factor in Eq. (34) to be constant, and thus physically irrelevant. If it is not feasible to control the density without inducing flows, then we must conclude that the density must be constant in time, and thus so |u| is constant in both space and time. This conclusion can also be drawn from the chemical potential µ of the BEC. If the BEC is stationary, then µ is constant. Since we can conclude that |u| is constant. Using Eq. (27), we can see that all of these restrictions on the density, speed of sound and flows are only achievable by balancing the external potential V and internal interaction strength λ to modulate the speed of sound while leaving the density constant in time. Specifically, we must have It is well known that the interaction strength in a BEC can be modulated with an external magnetic field around a Feshbach resonance (see for example [36]). Defining and for "unperturbed" interaction strength λ 0 and external potential V 0 , we must have From the definitions of c s and c 0 in Eqs. (22) and (23) it is straightforward to show that the interaction strength perturbation must have the form and thus where The density may still vary over space, with a shape determined by the "unperturbed" external potential V 0 as usual. Physically, this would result in a BEC cloud "not moving" in time (no flows, fixed density distribution), but with a carefully balanced trapping potential and applied magnetic field changing the speed of sound. This is somewhat analogous to modulating the refractive index in a dielectric, a scheme which has also been explored for its applications in analogue gravity (for example, in [7,37]). It should be noted that implementing the conditions presented in this Section does not necessarily result in an exact simulation, as the effective metrics in Eq. (34) and Eq. (28), with the above speed of sound perturbation, differ by a conformal factor, as in the simulation of various black hole geometries in [10,13,38]. Explicitly, in our case, The conformal factor will be discussed further in Section III C. While a GW is a coordinate independent physical effect, the simulation presented here reproduces elements of a metric in a particular coordinate system, and is thus not a coordinate independent solution.

Non-relativistic limit
In the explicitly non-relativistic limit, the spatial flows are much slower than the speed of light, so u 0 → c. The interaction strength must also be weak, so c 0 ≪ c. From the definition of c s in Eq. (22), it is clear that we must have c s ≪ c. In this regime, the phononic regime dispersion relation condition in Eq. (13) reduces to We also assume that the term 2 ∂ 2 t φ/mc 2 can be neglected, i.e. the excitation energy of each boson is much smaller than its mass energy. The equation governing the evolution of the field φ becomes the Gross-Pitaevskii equation, where the external potential and interaction strength are related to those defined in Section II B by In the non-relativistic limit, these are given by B. GW effect simulation The goal of this section is to simulate the acoustic metric given in [25], to test the metrological scheme proposed in [25,26]. This metric has the form where ρ 0 , c s0 and v 0 are the properties of the mean field φ "unperturbed" by a simulated GW. In [25], the BEC is considered to be at rest (v 0µ = −cδ 0 µ ) in the TT frame, so this is the condition that we will simulate with. Such a simulation can also be done with a BEC at rest in Fermi normal coordinates, which is not the initial condition considered in [25], and this solution is presented in the Appendix.

Acoustic metric with background GW
The metric perturbation h in Fermi normal coordinates is given in Eq. (8). Going from the TT frame to Fermi normal coordinates, the flow velocity transforms as Hence, the acoustic metric has the form This is the form of the acoustic metric that we will simulate.

Acoustic metric with simulated GW
To simulate the metric in Eq. (53), we perturb the bulk properties of the BEC in Fermi normal coordinates (i.e. lab frame). When constructing the simulation, we consider the background metric to be the flat Minkowski metric in Eq. (1), with no GW. Let the density, speed of sound and flow velocity be respectively described as and with normalisation so to first order, v 0 = c. With this flow velocity, the acoustic metric is to first order in ǫ. Comparison of Eq. (53) and Eq. (59) suggests that the velocity perturbation functions should take the form and the speed of sound perturbation must be As in Section III A, these conditions result in a conformal simulation, with the conformal factor G (sim) If all components of both effective metrics are to match, there is no way to avoid this conformal factor. The form of the density perturbation ρ 1 required to implements Eqs. (60) and (61) will be calculated in section III B 3. The conformal factor will be discussed further in section III C.

Bulk properties for simulation
To implement the normalised velocity profile given above, we must calculate the restrictions placed on the other bulk properties of the condensate. From Eq. (25), we can derive the velocity normalisation required for irrotational flow. We find that where x represents all spatial dimensions (x, y and z). As in Section III A, in the limit of ǫ → 0, the BEC is stationary. Hence, the chemical potential is constant, so |u| 0 is constant in both space and time. Similarly, from these results and the continuity equation, Eq. (B16), we can derive the form of the density and its perturbation. We find that the density and its perturbation have the form and where α (x) is some arbitrary function of integration, encoding the spatial shape of the BEC cloud. It should be noted that the results of this section are not fundamental restrictions on the bulk properties on the BEC; rather, they are conditions that must be imposed in an experiment to facilitate the implementation of the desired flow velocities and speed of sound.
As in Section III A, all of these conditions cannot be satisfied with an arbitrary interaction strength λ. Taking the same approach as above, we define a "perturbed" interaction strength In general from Eqs. (22) and (23), we must have With the results of this section, this expression becomes Using Eq. (27), this corresponds to an external potential

Static bulk solution
Consider a BEC trapped in a uniform box potential with infinite potential walls. In such a case, the density of the BEC is approximately constant in space everywhere inside the box, apart from a region close to the boundaries of the trap, where the density goes to zero. The width of this boundary region is given by the healing length defined above in Eq. (15). However, as stated in the motivation for the definition of Eq. (15), we are interested in perturbations whose wavelength far exceeds the healing length. Hence, for the perturbations we are considering, in a uniform box potential, we can assume constant density everywhere. This is also assumed in the detector proposal [25]. As in Section III A, it seems most reasonable to require that ρ 0 is constant in time, and thus |u| 0 is also. In this case, the perturbed bulk properties required to simulate a gravitational wave derived above can be simplified somewhat. Applying these conditions, we find implemented with Using the definitions and conditions presented in Section III A 3, the non-relativistic limit of these potential and interaction strength perturbations are where k µ is the wave-vector of the GW, andh µν is the trace reversed perturbation defined for metric perturbation h µν and Minkowski metric η µν as defined above. This follows simply from the Riemann tensor for a GW spacetime [32] R αµβν = 1 2 (h αν,µβ + h µβ,να − h µν,αβ − h αβ,µν ) .
In the TT gauge, the elements of the Weyl tensor have simple forms such as which can be measured by the detector and compared against experimental parameters of the simulation.

IV. EXAMPLES OF GW SOURCES
A. Non-axisymmetric neutron star Rotating neutron stars are one of the strongest predicted sources of continuous GWs [39]. Any imperfections in the symmetry of the mass distribution of a neutron star generate gravitational radiation as the star spins. The simplest case of a non-axisymmetric neutron star spinning down has strain functions of the form [40] ǫh where ι is the inclination of the neutron star's rotation axis to the line of sight, the phase evolution is for rotation frequency f /2 and reference time t 0 , and the amplitude h 0 is with ellipticity where I ii is the moment of inertia of the neutron star about some i axis, d is the distance to the neutron star and G is Newton's gravitational constant. This coordinate system is defined such that the axis of rotation is parallel to the z axis. On the time scale of a detection event, the frequency is constant to very good approximation, so terms in ∂ t f in the phase are ignored [40]. The signal emitted by such a neutron star can be directly simulated with the interaction and external potential perturbations B. Compact binary coalescence The first direct experimental proof of the existence of GWs was recently reported by the LIGO collaboration in [17], with the measurement of the GW signature of the final moments of a compact binary inspiral involving two black holes. These black holes were approximately 29 and 36 times the mass of the sun respectively, and 3 solar masses in energy was radiated in the form of GWs in the inspiral and collision. The form of the emitted gravitational radiation during the collision in the "strong gravity regime" must be calculated numerically, but the radiation emitted during the well separated inspiral phase, and the ringdown after coalescence, has well known solutions.

Inspiral
During the inspiral of a compact binary system, while the two compact objects are still well separated, the gravitational radiation far from the binary system has the form [32] ǫh where ι is the inclination axis of inspiral axis to detector, M = M 1 + M 2 and µ = M 1 M 2 /M 2 for the two masses M 1 and M 2 , d is the distance from the inspiral barycentre to the detector, and with some reference time t 0 . This is a sinusoidal signal whose amplitude and frequency increase as the time t reaches the reference time t 0 , i.e. the time of collision. This is the characteristic "chirp" observed by the LIGO collaboration in [17][18][19]. To directly simulate this metric, the corresponding interaction and external potential perturbations

Ringdown
After a binary system with sufficient mass to form a black hole has collided and coalesced, the resulting black hole rotates due to conservation of angular momentum. The ringdown of the coalesced object into a stable rotating black hole can thus be modelled as a perturbed Kerr black hole. The simplest single-mode ringdown of a Kerr black hole has strain functions of the form [41] ǫh where as above, ι is still the inclination angle of rotation axis to the detector, d is the distance from the source to the detector, Q is the "quality factor" fitted numerically with for spin parameterâ = cS/GM 2 , with spin angular momentum S. The GW amplitude A is given by where F (Q) = 1 + 1/4Q 2 , g (â) = 1.5251 − 1.1568 (1 −â) 0.1292 and ε is the fraction of the black hole mass radiated away. Functionally, this is a decaying sinusoid of constant frequency.
The corresponding interaction and external potential perturbations for simulation are

V. REDUCTION TO 1+1
In this Section, we restrict ourselves to an effective 1-dimensional field to compare to earlier work in [16]. In an effective 1 + 1 dimensional spacetime, the GW metric reduces to To simulate this, the speed of sound is chosen as and a time scaling of results in a simulation In 1 + 1 dimensions, the equations of motion are conformally independent, so this is an exact simulation. Following the same procedure as in Section III A, we require that the flow velocity normalisation is completely constant and conclude that the density is constant in time. The interaction strength and external potential perturbations required to implement this are then We must stress that this is an effective 1 + 1 dimensional theory, and care must be taken when dealing with the actual field dynamics. Although this seems to work at the level of the metric, a naive suppression of the remaining spatial dimensions cannot be done due to the fact that the conformal factor is dimensionally dependent, and diverges when the number of spatial dimensions is exactly 1 [42]. Nevertheless, as long as the system is sufficiently constrained in the extra dimensions, e.g. in a highly elongated trap, a well-behaved effective 1 + 1 dimensional system can always be constructed.

VI. CONCLUSION
We have shown how to simulate a GW spacetime in 3 + 1 dimensions for quantum excitations of a BEC, up to a conformal factor, as well as simulating the acoustic metric used in [25] to propose a GW detector. By making use of the "gauge freedom" of the GW metric corresponding to a linearised coordinate transformation, we chose a frame in which the metric perturbation could be simulated by perturbing the speed of sound in the BEC. We then examined the restrictions this places on other bulk properties through the continuity equation and experimental limitations, and calculated the external and interaction potential perturbations needed to implement such a simulation in the lab. Although the simulated metric is related to the target metric by a non-constant conformal factor, we show that there are still useful properties that can be measured and tested in an experiment. We also give explicit expressions for the simulation of GWs from various sources. This work generalises the results of [16] and presents a complementary approach to simulation in effectively 1 + 1 dimensional BECs.
The results presented here can also be derived in the context of an explicitly nonrelativistic treatment of a BEC, such as that derived in [42]. In a non-relativistic BEC, phonons on the BEC still propagate on a Lorentzian effective spacetime described by an acoustic metric, but this metric is necessarily spatially conformally flat. We consider a BEC in a covariant formalism in this paper to match the approach of [25][26][27] for the simulation in Section III B, and to express the interaction of a BEC with GWs in a natural way. As explained in Section II B, the simulation of GWs presented in this paper do not rely on the relativistic nature of the BEC or any relativistic effects, nor do the perturbations to the external and interaction potentials disappear in the non-relativistic limit.
We have studied GWs in the context of perturbations around a flat spacetime metric and assuming GWs to be far outside the source. Other interesting simulations could involve GWs propagating on curved backgrounds, such as black holes [9,10] or during inflation [43], or in strong-field regimes. Furthermore, since phonons are quantum quasi-particles, this opens up the possibility of studying predicted effects of quantum field theory in curved spacetime, such as how a GW may affect the entanglement of quantum systems, a phenomena that is utilised in the BEC GW detector proposed in [25]. This, therefore, also presents a potential, and fully configurable, testing environment for this GW detector metrological scheme. To obtain a full simulation of the GW detector, we need a better understanding of the effect of the GW on the bulk of the BEC. As mentioned in the conclusion of [16], an experimental simulation of the effect of a large amplitude GW and subsequent detection of phonons would also be a proof-of-concept demonstration of the generation of phonons by GWs as predicted in [25]. As explained in Section III B, it is possible to simulate the effect of a GW starting with no flows in Fermi normal coordinates. This solution is presented here. As above, let the density, speed of sound and flow velocity normalisation be respectively described as and where ǫ is the small parameter defined above. ρ 0 , c s0 , and |u| 0 are bulk properties of the BEC in the absence of a simulated GW.

Effective metric with background gravitational wave
In Fermi normal coordinates, we consider the flows in the BEC to be v µ = −v 0 δ 0 µ . With the normalisation equation (Eq. (24)), we can determine the function v 0 as Then, with a background GW, the acoustic metric has the form To simulate the effect of a GW, we perturb the bulk properties of the BEC. The acoustic metric is

Simulation
Comparison of Eq. (A6) and Eq. (A8) suggests that, to simulate a background GW in these coordinates, the speed of sound should be modulated as As with the simulation presented in the main body of the paper, this simulation differs from an exact simulation by a conformal factor;

a. Bulk properties for simulation
To implement the normalised velocity profile given above, we must calculate the restrictions placed on the other bulk properties of the condensate. From Eq. (25), and from Eq. (B16), As in Section III A, it seems most reasonable to require that ∂ t ρ 0 = 0 and so |u| 0 is completely constant. Defining a "perturbed" interaction strength as the results of this Section and Eqs. (22) and (23) are simultaneously satisfied if This can be implemented together with the external potential Appendix B: Acoustic metric with general background metric In [13,34], the acoustic metric is derived for a flat Minkowski background metric. Since we require the same for a general background metric, in this appendix we extend the acoustic metric to the case where the background metric is not necessarily flat.

a. Lagrangian
The Lagrangian density for an interacting massive complex scalar fieldΦ on a (in general curved) background with metric g µν may be written as where m is the mass, the external potential V is generally a function of space and time, and the interaction potential U depends on coupling constants λ i which are also in principle functions of space and time. The background metric g µν cannot be completely general; we restrict ourselves to spacetimes with sufficiently weak curvature such that Bose-Einstein condensation can still be well defined. Further restrictions on the metric will be given in Section B 2 b. The interaction potential U can be expanded as We will consider only the first term of U corresponding to two-particle interactions, and ignore further terms corresponding to three or more particle interactions. For notational convenience, we will drop the label on λ 2 so The Euler-Lagrange equation forΦ † is g − and g is the determinant of g µν .

b. Approximations
We now let this fieldΦ represent a Bose-Einstein condensate and make the Bogoliubov approximation to separate the "condensed fraction" of the field φ from a small "uncondensed fraction"ψ. This is done multiplicatively aŝ to simplify the equation forψ later. As part of the Bogoliubov approximation, we say that where · is a non-equilibrium average. Taking the average of Eq. (B4), We now take the Popov approximation and require that the density of excited atoms is much smaller than the density of mean-field atoms, i.e. ψ †ψ ≪ 1.
This results in a non-linear Klein-Gordon-like equation for the mean field φ: This is a curved space-time generalisation of the Gross-Pitaevskii equation. In flat space-time (where the metric is the Minkowski metric η µν defined in Eq. (1)) and in the non-relativistic limit, we can replace φ with a lower energy field φ = ϕe imc 2 t/ (B12) and take the limit of c → ∞. Assuming that the energy of excitations in ϕ is sufficiently low such that we can ignore terms of order ∂ 2 t ϕ, the remaining terms of Eq. (B11) have the form which is the usual time dependent Gross-Pitaevskii equation.

c. Continuity and velocity normalisation equations
If the mean field φ is written in the Madelung representation φ = √ ρe iθ and defining a flow velocity where andT is a generalised kinetic operator, which reduces to the standard kinetic energy operator T = − ( 2 /2m) ∇ 2 for constant ρ in the non-relativistic flat space-time limit. Note that we require the solution to Eq. (B11) to solve Eq. (B18) but not vice-versa, as we are neglecting the back-reaction ofψ on φ. Taking the equivalent equation to Eq. (B18) forψ † and combining these to eliminateψ † , we find i u µ ∂ µ +T ρ 1 c 2 0 −i u µ ∂ µ +T ρ + 2mT ρ ψ = 0.
It is important to note that although Eq. (B18) implies Eq. (B21), the converse is not true.

b. Relative term strength
The phonon equation Eq. (B21) can be expanded into four terms aŝ andT We make an eikonal approximation, where and the corresponding relations for variations in space as in the flat space case, but also with the corresponding relations for variations in space. Note that Eq. (B28) restricts the curvature of the metric with respect to the phonon mode frequencies. For linearised gravity and realistic phonon frequencies, this will always hold. Additionally, following [13,35] we consider small momenta in the phononic regime, such that the dispersion relation is linear and terms quartic in k can be neglected. With these approximations,T 2 andT 3 are negligible in comparison toT 1 andT 4 , so we are left with