Exploring the unification of quantum theory and general relativity with a Bose–Einstein condensate

1.1. Motivation and background

1.2. Experimental approaches

1.3. Outline

Here we consider how GQSR can arise due to a conflict between the superposition principle of QT and equivalence principle of GR⁹ . More detail can be found in [4, 6, 34]. Also, see [3, 34] for how the same proposed state reduction can be derived using the principle of covariance rather than the principle of equivalence.

Let us first consider a simple situation of a tabletop quantum experiment where it is required that the Earth's gravitational field is to be taken into consideration (see figure 1). There are basically two different ways to incorporate the Earth's field in this experiment (which is to be considered as constant, both spatially and temporally, and to be treated non-relativistically). The first, the Newtonian perspective, would simply be to incorporate a term in the Hamiltonian representing the Newtonian potential (this being the normal prescription that most physicists would adopt), and use standard Cartesian coordinates (x, y, z, t), or rather $({\boldsymbol{r}},t)$ in three-vector form. The second, the Einsteinian perspective, would be to adopt a freely falling reference system $({\boldsymbol{R}},T)$, in accordance with which the Earth's gravitational field vanishes. The relation between the two is:

$$\begin{eqnarray}&&{\boldsymbol{R}}={\boldsymbol{r}}-\displaystyle \frac{1}{2}{t}^{2}{\boldsymbol{a}},\hspace{1cm}T=t,\end{eqnarray} \tag{ 2 }$$

where the constant three-vector ${\boldsymbol{a}}$ is the acceleration due to the Earth's gravity. We denote the wavefunction in the $({\boldsymbol{r}},t)$ system, using the Newtonian perspective, by ψ, whereas for the $({\boldsymbol{R}},T)$ system, using the Einsteinian perspective, we use Ψ. For a free particle of mass m, we have, according to the Newtonian perspective, the Schrödinger equation

$$\begin{eqnarray}&&{\rm{i}}{\hslash }\displaystyle \frac{\partial \psi }{\partial t}=-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}^{2}\psi -m{\boldsymbol{r}}.{\boldsymbol{a}}\psi ,\end{eqnarray} \tag{ 3 }$$

whereas, according to the Einsteinian perspective

$$\begin{eqnarray}&&{\rm{i}}{\hslash }\displaystyle \frac{\partial {\rm{\Psi }}}{\partial t}=-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}^{2}{\rm{\Psi }},\end{eqnarray} \tag{ 4 }$$

the operator ∇² being the same in both coordinate systems. To get consistency between the two perspectives, we need to relate ψ to Ψ by a phase factor [4–6, 45–47]

$$\begin{eqnarray}&&{\rm{\Psi }}={{\rm{e}}}^{{\rm{i}}\tfrac{m}{{\hslash }}(\tfrac{1}{6}{t}^{3}{a}^{2}-t{\boldsymbol{r}}.{\boldsymbol{a}})}\psi .\end{eqnarray} \tag{ 5 }$$

Zoom In Zoom Out Reset image size

For a quantum experiment involving many particles of total mass $\bar{m}$ and center of mass $\bar{{\boldsymbol{r}}}$ (or $\bar{{\boldsymbol{R}}}$ in the Einstein system), this generalizes to

$$\begin{eqnarray}&&{\rm{\Psi }}={{\rm{e}}}^{{\rm{i}}\tfrac{\bar{m}}{{\hslash }}(\displaystyle \frac{1}{6}{t}^{3}{a}^{2}-t\bar{{\boldsymbol{r}}}.{\boldsymbol{a}})}\psi .\end{eqnarray} \tag{ 6 }$$

Since the difference between the Newtonian and Einsteinian perspectives is merely a phase factor, one might form the opinion that it makes no difference which perspective is adopted. Indeed, the famous experiment by Colella, Overhauser and Werner [48] (see also [49–51]) performed originally in 1975, did provide some impressive confirmation of the consistency (up to a certain point) of QT with Einstein's principle of equivalence.

However, it is important to note that the phase factor that is encountered here is not at all harmless, as it contains the time-dependence involved in the term

$$\begin{eqnarray}&&\displaystyle \frac{1}{6}{t}^{3}{a}^{2},\end{eqnarray} \tag{ 7 }$$

in the exponent, which affects the splitting of field amplitudes into positive and negative frequencies. In other words, the Einsteinian and Newtonian wavefunctions belong to different Hilbert spaces, corresponding to different quantum field theoretic vacua. In fact, this situation is basically just a limiting case of the familiar relativistic (Fulling-Davies-)Unruh effect [45, 46, 52–55], where in a uniformly accelerating (Rindler) reference frame, we get a non-trivial thermal vacuum of temperature

$$\begin{eqnarray}&&\displaystyle \frac{{\hslash }a}{2\pi {kc}},\end{eqnarray} \tag{ 8 }$$

where a is the magnitude of acceleration, k being Boltzmann's constant and c, the speed of light. In the current situation, we are considering the Newtonian approximation $c\to \infty$, so the temperature goes to zero. Nevertheless, as a direct calculation shows, the Unruh vacuum actually goes over to the Einsteinian one in the limit $c\to \infty$, in agreement with what has been shown above, and is thus still different from the Newtonian one even though the temperature difference goes to zero in this limit.

At this stage, we could still argue that it makes no difference whether the Newtonian or Einsteinian perspective is adopted, so long as one sticks consistently to one perspective or the other overall (since the formalism is maintained within a single Hilbert space). However, the situation becomes radically different when one considers the gravitating body, in this example the Earth, to be in a quantum superposition between pairs of states in which the gravitational fields differ. If we were to adopt the Newtonian perspective for our quantum experiment then we would encounter no problem with the formalism of QT, the standard linear framework of unitary evolution applying as well to the Newtonian gravitational field as it does to electromagnetism, or to any other standard field theory of physics. But it is another matter altogether if we insist on adopting the Einsteinian perspective. Our standpoint here is that, owing to the enormous support that observations have now provided for GR in macroscopic physics, one must try to respect the Einsteinian perspective as far as possible, in quantum systems, especially in view of the foundational role that the principle of equivalence has for GR (see [4–6]).

Let us now replace the Earth with a small rock and try to imagine the quantum description of the physics taking place in some small region in the neighborhood of the rock, where we consider that the rock can persist for some while in a superposition of two different locations, and we label the respective states as $| L\rangle$ and $| R\rangle$. We are not now trying to compare the Einsteinian perspective with a Newtonian one, since our point of view will be that the latter is not relevant to our current purposes, as we regard the Einsteinian perspective to be closer to nature's ways. Instead, we attempt to adopt an Einsteinian perspective for a quantum experiment in the neighborhood of the rock that is in a quantum superposition of two locations, $\alpha | L\rangle +\beta | R\rangle$. What we now have to contend with is a superposition of two different Einsteinian wavefunctions for the quantum experiment, each inhabiting a Hilbert space that will turn out to be incompatible with the other.

However, the preceding discussion does not hold exactly, because for each of the two components of the superposition of rock locations $| L\rangle$ and $| R\rangle$, the gravitational field of the rock is not completely uniform. Nevertheless, we shall consider, first, that we are examining the nature of the wavefunction in some spatial region that is small by comparison with the rock itself, so that we can assume that the gravitational field of each component of the superposition can be taken to be spatially uniform to a good approximation. Adopting the Einsteinian perspective, what we are now confronted with is the fact that the gravitational acceleration fields for the two rock locations will be different from each other, so that the difference between these local acceleration fields ${\boldsymbol{a}}$ and ${\boldsymbol{a}}^{\prime}$ will lead to a difference between the Einsteinian vacuum for each rock location. In the neighborhood of each spatial point, there will be a phase difference between the two states of our quantum experiment that are in superposition:

$$\begin{eqnarray}&&{{\rm{e}}}^{{\rm{i}}\tfrac{\bar{m}}{{\hslash }}\left(\displaystyle \frac{1}{6}{t}^{3}{\left({\boldsymbol{a}}-{\boldsymbol{a}}^{\prime} \right)}^{2}-t\bar{{\boldsymbol{r}}}.({\boldsymbol{a}}-{\boldsymbol{a}}^{\prime} )\right)}.\end{eqnarray} \tag{ 9 }$$

Although the presence of the $\tfrac{1}{6}{t}^{3}{\left({\boldsymbol{a}}-{\boldsymbol{a}}^{\prime} \right)}^{2}$ term tells us, strictly speaking, that when ${\boldsymbol{a}}\ne {\boldsymbol{a}}^{\prime}$, the superposition is illegal (the states belonging to different Hilbert spaces), we adopt the view that this incompatibility takes some time to cause trouble (as would eventually become manifest in divergent scalar products, etc). The idea is that in order to resolve this incompatibility of Hilbert spaces, the superposed state must eventually reduce to one alternative or the other, this incompatibility taking a while to build up. We compare the troublesome term $\tfrac{1}{6}{t}^{3}{\left({\boldsymbol{a}}-{\boldsymbol{a}}^{\prime} \right)}^{2}$ with the harmless one $t\bar{{\boldsymbol{r}}}.({\boldsymbol{a}}-{\boldsymbol{a}}^{\prime} )$, the latter ($\times \bar{m}/{\hslash }$) being linear in t and therefore not altering the vacuum but, in effect, just corresponds to incorporating the Newtonian gravitational potential term into the Hamiltonian. We take the view that so long as t is small enough, the trouble arising from t³ remains insignificantly small, where the measure of this smallness comes from comparing $\tfrac{1}{6}{t}^{3}{\left({\boldsymbol{a}}-{\boldsymbol{a}}^{\prime} \right)}^{2}$ with the harmless $t\bar{{\boldsymbol{r}}}.({\boldsymbol{a}}-{\boldsymbol{a}}^{\prime} )$.Thus, we take the coefficient $\tfrac{1}{6}{t}^{3}{\left({\boldsymbol{a}}-{\boldsymbol{a}}^{\prime} \right)}^{2}$ as some kind of measure of the propensity for the state to reduce, as a contribution to the overall reduction process. To get our measure of total error, or 'uncertainty' Δ, we integrate this expression over the whole of (coordinate) 3-space:

$$\begin{eqnarray}&&{\rm{\Delta }}:= \gamma \int {\left({\boldsymbol{a}}-{\boldsymbol{a}}^{\prime} \right)}^{2}{{\rm{d}}}^{3}{\boldsymbol{r}},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&=\,\gamma \int {\left({\boldsymbol{\nabla }}\phi -{\boldsymbol{\nabla }}\phi ^{\prime} \right)}^{2}{{\rm{d}}}^{3}{\boldsymbol{r}},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&=\,\gamma \int \left[{\boldsymbol{\nabla }}(\phi -\phi ^{\prime} )\right].\left[{\boldsymbol{\nabla }}\phi -{\boldsymbol{\nabla }}\phi ^{\prime} \right]{{\rm{d}}}^{3}{\boldsymbol{r}},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&=\,-\gamma \int (\phi -\phi ^{\prime} )\left({{\rm{\nabla }}}^{2}\phi -{{\rm{\nabla }}}^{2}\phi ^{\prime} \right){{\rm{d}}}^{3}{\boldsymbol{r}},\end{eqnarray} \tag{ 13 }$$

(assuming appropriate falloff at spatial infinity), where γ is some positive constant, and ϕ and $\phi ^{\prime}$ are the respective gravitational potentials for the states of the rock, where we are adopting a Newtonian approximation for estimating the required error (${\boldsymbol{a}}=-{\boldsymbol{\nabla }}\phi$ and ${\boldsymbol{a}}^{\prime} =-{\boldsymbol{\nabla }}\phi ^{\prime}$). By Poisson's formula (G being Newton's gravitational constant)

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}\phi =4\pi G\rho ,\end{eqnarray} \tag{ 14 }$$

we get

$$\begin{eqnarray}&&{\rm{\Delta }}=-4\pi G\gamma \int (\phi -\phi ^{\prime} )(\rho -\rho ^{\prime} ){{\rm{d}}}^{3}{\boldsymbol{r}},\end{eqnarray} \tag{ 15 }$$

where ρ and $\rho ^{\prime}$ are the respective mass densities of the two states, and we shall take these mass densities in the sense of expectation values for the respective quantum states. Using the formula

$$\begin{eqnarray}&&\phi ({\boldsymbol{r}})\equiv -G\int \displaystyle \frac{\rho ({\boldsymbol{r}}^{\prime} )}{| {\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} | }{{\rm{d}}}^{3}{\boldsymbol{r}}^{\prime} ,\end{eqnarray} \tag{ 16 }$$

we obtain [3]:

$$\begin{eqnarray}&&{\rm{\Delta }}=4\pi {G}^{2}\gamma \int \int \displaystyle \frac{\left[\rho ({\boldsymbol{r}})-\rho ^{\prime} ({\boldsymbol{r}})\right]\left[\rho ({\boldsymbol{r}}^{\prime} )-\rho ^{\prime} ({\boldsymbol{r}}^{\prime} )\right]}{| {\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} | }{{\rm{d}}}^{3}{\boldsymbol{r}}{{\rm{d}}}^{3}{\boldsymbol{r}}^{\prime} .\end{eqnarray} \tag{ 17 }$$

Defining E_G := Δ/G, we have a quantity that is proportional to (depending on the value of γ) the gravitational self-energy of the difference between the mass distributions of each of the two states

$$\begin{eqnarray}&&{E}_{G}=4\pi \gamma \int (\phi -\phi ^{\prime} )(\rho ^{\prime} -\rho ){{\rm{d}}}^{3}{\boldsymbol{r}}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&=4\pi G\gamma \int \int \displaystyle \frac{\left[\rho ({\boldsymbol{r}})-\rho ^{\prime} ({\boldsymbol{r}})\right]\left[\rho ({\boldsymbol{r}}^{\prime} )-\rho ^{\prime} ({\boldsymbol{r}}^{\prime} )\right]}{| {\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} | }{{\rm{d}}}^{3}{\boldsymbol{r}}{{\rm{d}}}^{3}{\boldsymbol{r}}^{\prime} .\end{eqnarray} \tag{ 19 }$$

The quantity Δ can be considered as a measure of a limitation to regarding the quantum superposition of the rock $\alpha | L\rangle +\beta | R\rangle$ as being a stationary state, in accordance with principles of GR. Thus, we may take it to be a reasonable inference from general-relativistic principles to regard Δ⁻¹ as providing some kind of measure of a limit to the length of time that the superposition might persist, the shorter that this time-scale should presumably be, the larger the value Δ is found to have. This conclusion comes from considerations of GR, as applied simply to the notion of a quantum superposition of states, no consideration of quantum dynamics being involved except for the quantum notion of stationarity. Moreover, no actual measure of a time-scale for a 'lifetime' of the superposition has yet been provided by these considerations.

However, a significant clue is provided by Heisenberg's time-energy uncertainty principle, where we note that the quantity E_G = Δ/G is an energy. In QT, the lifetime of an unstable particle is reciprocally related to an energy uncertainty, which can be regarded as a manifestation of Heisenberg's time-energy uncertainty principle. In a similar way, we propose that E_G should be treated as a fundamental uncertainty in the energy of the superposition $\alpha | L\rangle +\beta | R\rangle$. We then take the view that the 'energy uncertainty' E_G is reciprocally related to a lifetime of this superposition between the states $| L\rangle$ and $| R\rangle$, and we can, therefore, regard the macroscopic superposition as having an average lifetime τ that is roughly given by

$$\begin{eqnarray}&&\tau \sim \displaystyle \frac{{\hslash }}{{E}_{G}},\end{eqnarray} \tag{ 20 }$$

upon which time (on average) the state $\alpha | L\rangle +\beta | R\rangle$ spontaneously 'decays' into one or the other of $| L\rangle$ or $| R\rangle$. This decay process cannot be derived from considerations of QT alone, and instead we are assuming the invocation of a higher theory from which QT and GR are limiting cases. The energy uncertainty in (20) arises due to a conflict between the general-relativistic and quantum principles that are being appealed to in relation to the description of stationary gravitational fields in quantum linear superposition, and there would be no need for such an energy uncertainty if we had just assumed a Newtonian description of gravity without the philosophy of GR. Similarly, if we had considered a contribution from the electromagnetic interaction of a (say charged) rock in addition to its gravitational field, then there would be no conflict with QT from electromagnetic effects (there being no equivalence principle for electromagnetism) and we would not be led to consider any energy uncertainty from electromagnetic effects contributing to the decay of the state to either $| L\rangle$ or $| R\rangle$.

Taking the analogy with particle decay further, we could assume that the probability of, a presumed spontaneous, state reduction is an exponential function of time t:

$$\begin{eqnarray}&&{P}_{s}(t)={{\rm{e}}}^{-t/\tau }={{\rm{e}}}^{-{E}_{G}t/{\hslash }},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{P}_{d}(t)=1-{{\rm{e}}}^{-t/\tau }=1-{{\rm{e}}}^{-{E}_{G}t/{\hslash }},\end{eqnarray} \tag{ 22 }$$

where P_s(t) and P_d(t) are, respectively, the probability of survival and decay of the superposition state. Here we are assuming, as with particle decay, that the decay is memoryless, which would seem the simplest assumption for describing the decay process given that there is not, at present, a full theory. Equation (22) illustrates that it should not be necessary to wait for a time τ = E_G/ℏ in order to observe collapse of the wavefunction, and we can estimate how often collapse will occur at a given time t without having to appeal to a full dynamical model.

A few clarifying remarks should be made on our above derivation of E_G. We have considered a rock to be in a superposition of two locations (similar to Schrödinger's cat being in two locations). However, rather than a rock (which we have assumed to be a continuous mass distribution), we could have considered a point-like object in the superposition of two locations, which, naively, would be a superposition of delta-function wavefunctions in position space. This would lead to an infinite value for E_G. The problem is that a delta function is not a stationary solution to the Schrödinger equation since the position wavefunction would instantly spread out (there is infinite uncertainty in momentum). However, the stationary solution of the Schrödinger equation would be that where the state is spread out over the Universe (there is infinite uncertainty in position), which is clearly not satisfactory either. One might imagine that, in a full theory of GQSR, a spreading state keeps reducing by GQSR. For now, a systematic procedure would be to modify the Schrödinger equation and use the Schrödinger–Newton equation [56] to obtain stationary states [6, 57]. For a point-like object, the stationary solution is then a 'smeared-out' delta function, and the position is no longer defined at a point. To calculate E_G for this state, we could think of each 'smeared-out' delta function as representing point-like objects in a superposition of continuous positions with differing weights, and follow the procedure outlined with equations (10)–(19). In the continuous limit, this results in the same expression for E_G as before, (19), with $\rho ({\boldsymbol{r}})$ and $\rho ^{\prime} ({\boldsymbol{r}})$ representing the average mass density of the defined stationary states, which would be $m| \psi ({\boldsymbol{r}}){| }^{2}$ and $m| \psi ^{\prime} ({\boldsymbol{r}}){| }^{2}$ in our case, where ψ and $\psi ^{\prime}$ are the normalized wavefunctions of the stationary states, and m is the mass of the object. The superposition state should then decay into one of these stationary states as outlined in (20)–(22). For our rock, we have assumed that the stationary states $| L\rangle$ and $| R\rangle$ should be close to the 'classical' rock states since we do not see rocks spreading out across the Universe¹⁰ . The mass profile $\rho ({\boldsymbol{r}})$ of the stationary state of the rock should then be close to its classical mass distribution, which we have approximated as a continuous mass distribution.

The above calculation of E_G, (10)–(19), has been carried out entirely within the framework of Newtonian mechanics since we are considering the masses involved as being rather small and moving slowly, so that general-relativistic corrections can be ignored to a very good approximation. We can, therefore, also just use Schrödinger's equation rather than, for example, the full framework of quantum field theory. However, the principles of GR still apply to gravity, such as the equivalence principle, and the fact that E_G is to be regarded as an energy uncertainty is coming from considerations of general-relativistic principles and QT. The use of Newtonian mechanics for calculating an expression for E_G, while nevertheless retaining much of the basic philosophy of Einstein's theory, is perhaps most clearly expressed with the Newton–Cartan theory of gravity [4].

2.1. Gravitationally-induced state reduction for uniform spherical mass distributions

In order to get an impression of the role of E_G, we can first think of the case of a solid spherical ball of radius R, made from some uniform massive material of total mass M, where the ball is in a superposition of two locations, differing by a distance b. The quantity E_G in this case is (see, for example [6], and appendix A):

$$\begin{eqnarray}{E}_{G}=\left\{\begin{array}{cl}\displaystyle \frac{6{{GM}}^{2}}{5R}\left(\displaystyle \frac{5}{3}{\lambda }^{2}-\displaystyle \frac{5}{4}{\lambda }^{3}+\displaystyle \frac{1}{6}{\lambda }^{5}\right) & \mathrm{if}\ 0\leqslant \lambda \leqslant 1,\\ \displaystyle \frac{6{{GM}}^{2}}{5R}\left(1-\displaystyle \frac{5}{12\lambda }\right) & \mathrm{if}\ \lambda \geqslant 1,\end{array}\right.\end{eqnarray} \tag{ 23 }$$

where λ := b/(2R), and we have taken γ = 1/(8π) in (18). See figure 2(a) for an illustration of E_G, and figure 2(b) for the rate of change of E_G, with separation b. The only point of particular relevance is the fact that, for a displacement such that the two spheres are touching, the value of E_G is already nearly $\tfrac{2}{3}$ of the value it would reach for a displacement all the way out to infinity. Thus, for a uniformly solid body like this, we do not gain much by considering displacements in which the two instances of the body are moved apart by a distance of more than roughly its own diameter.

Zoom In Zoom Out Reset image size

2.2. Gravitationally-induced state reduction for uniform spheroidal mass distributions

The above case of two uniform spherical mass distributions is that which is generically considered in the literature, apart from a study of two uniform cubes at only very small displacement b [32]¹¹ . In this section, we generalize to uniform spheroidal mass distributions and consider whether this can lead to an increase in E_G, and thus a faster rate of state reduction.

Now that we no longer have spherical symmetry, there are various configurations for the displacement of the spheroids. Here we consider four possible configurations: (a) an oblate spheroid displaced along its symmetry axis, (b) a prolate spheroid also displaced along its symmetry axis, (c) an oblate spheroid displaced along an equatorial (semi-major) axis, and (d) a prolate spheroid also displaced along an equatorial (now semi-minor) axis. See figure 3 for a visual illustration of all these configurations. Although analytical solutions can be obtained for a general expression of E_G for these cases (i.e. E_G for a general equatorial or polar displacement b between the spheroid states), the results are rather cumbersome and here we instead provide the results for the two cases (a) and (b) in the limit of high ellipticity e (see appendix D for more detail). Defining as $\epsilon := \sqrt{1-{e}^{2}}$, such that:

$$\begin{eqnarray}\epsilon := \left\{\begin{array}{cl}c/a & \mathrm{for}\ {\rm{an}}\,{\rm{oblate}}\,{\rm{spheroid}}\,(a\gt c),\\ a/c & \mathrm{for}\ {\rm{a}}\,{\rm{prolate}}\,{\rm{spheroid}}\,(c\gt a),\end{array}\right.\end{eqnarray} \tag{ 24 }$$

where a and c are the equatorial and polar radii respectively¹² , then, when e ≈ 1 we have $\epsilon \ll 1$. For the extreme ($\epsilon \ll 1$) prolate (spindle-like) spheroid in configuration (b), we find that, to first order in :

$$\begin{eqnarray}{E}_{G}=\left\{\begin{array}{cl}\displaystyle \frac{6{{GM}}^{2}}{5c}\left(A/4-B\mathrm{ln}\epsilon \right) & \mathrm{if}\ 0\leqslant b\leqslant 2c,\\ \displaystyle \frac{6{{GM}}^{2}}{5c}\left(C/4-\mathrm{ln}\epsilon \right) & \mathrm{if}\ b\geqslant 2c,\end{array}\right.\end{eqnarray} \tag{ 25 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}A & := & 2\lambda \left[1+\mathrm{ln}(1+\lambda )\right]+{\lambda }^{2}(20-21\mathrm{ln}2)-{\lambda }^{3}\left[21-10\mathrm{ln}\displaystyle \frac{{\left(1+\lambda \right)}^{2}}{4\lambda }\right]\\ & & +77{\lambda }^{4}-2{\lambda }^{5}\left[35+\mathrm{ln}\displaystyle \frac{1+\lambda }{4\lambda }\right]+20{\lambda }^{6},\end{array}\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&B:= 5{\lambda }^{2}-5{\lambda }^{3}+{\lambda }^{4},\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&C:= 4\mathrm{ln}2-11\lambda -4(1-5{\lambda }^{2}){\coth }^{-1}\lambda +2{\lambda }^{3}\left[({\lambda }^{2}-5)\mathrm{ln}\displaystyle \frac{{\lambda }^{2}}{1-{\lambda }^{2}}-1\right]\end{eqnarray} \tag{ 28 }$$

with λ now defined as λ := b/(2c).

Zoom In Zoom Out Reset image size

On the other hand, for the extreme ($\epsilon \ll 1$) oblate (pancake-like) spheroid in configuration (a), we find that E_G can be approximated by:

$$\begin{eqnarray}{E}_{G}=\left\{\begin{array}{cl}\displaystyle \frac{6{{GM}}^{2}}{5a}{ \mathcal A } & \mathrm{if}\ 0\leqslant b\leqslant 2c,\\ \displaystyle \frac{6{{GM}}^{2}}{5a}\left({ \mathcal C }/4-\epsilon \right) & \mathrm{if}\ b\geqslant 2c,\end{array}\right.\end{eqnarray} \tag{ 29 }$$

where A and C are defined as:

$$\begin{eqnarray}&&{ \mathcal A }:= 5{\beta }^{2}\left(\displaystyle \frac{1}{\epsilon }+\displaystyle \frac{\pi }{2}\right)-\displaystyle \frac{5{\beta }^{3}}{\epsilon }\left(\displaystyle \frac{\pi }{32}-\displaystyle \frac{1}{2\epsilon }\right)+\displaystyle \frac{{\beta }^{5}}{4{\epsilon }^{3}}\left(\displaystyle \frac{9\pi }{8}+\displaystyle \frac{1}{\epsilon }\right),\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{ \mathcal C }:= 2\pi +11\beta -2{\beta }^{3}-4(1+5{\beta }^{2}){\cot }^{-1}\beta +2{\beta }^{3}(5+{\beta }^{2})\mathrm{ln}\displaystyle \frac{1+{\beta }^{2}}{{\beta }^{2}},\end{eqnarray} \tag{ 31 }$$

with β := b/(2a). In figures 4(a) and (b), we provide E_G for the sphere and the above four spheroidal cases (a)–(d) when  = 0.5 (e = 0.87) and  = 0.01 (e = 0.999 95)¹³ . In all cases we take the volume and mass (and so density) of the objects to be the same. These figures illustrate that, for configurations (a) and (d), E_G can be greater than that of the sphere (with the same volume and density) at certain displacements, although, at infinite displacement the sphere has the greatest E_G. Indeed, at $b=\infty$, we find

$$\begin{eqnarray}&&{E}_{G}^{{\rm{sphere}}}=\displaystyle \frac{6{{GM}}^{2}}{5R},\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&{E}_{G}^{{\rm{prolate}}}=\displaystyle \frac{6{{GM}}^{2}}{5l}{\tanh }^{-1}e\equiv \displaystyle \frac{6{{GM}}^{2}}{5l}{\cosh }^{-1}\epsilon ,\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{E}_{G}^{{\rm{oblate}}}=\displaystyle \frac{6{{GM}}^{2}}{5l}{\sin }^{-1}e\equiv \displaystyle \frac{6{{GM}}^{2}}{5l}{\sec }^{-1}\epsilon ,\end{eqnarray} \tag{ 34 }$$

for the three cases of sphere, prolate and oblate, irrespective of how they are displaced with respect to each other, and l is the focal distance of the spheroids, which is $\sqrt{{a}^{2}-{c}^{2}}$ for the oblate, and $\sqrt{{c}^{2}-{a}^{2}}$ for the prolate spheroid. Equations (32)–(34) are valid for any value of e and between 0 and 1, and no constraints are placed on the size of volume and density. However, taking all the objects to have the same volume and mass, and assuming low ellipticity $e\ll 1$, we find

$$\begin{eqnarray}&&{E}_{G}^{{\rm{prolate}}}\approx \displaystyle \frac{6{{GM}}^{2}}{5R}\left(1-\displaystyle \frac{1}{45}{e}^{4}-\displaystyle \frac{64}{2835}{e}^{6}\,-\,\cdots \right),\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{E}_{G}^{{\rm{oblate}}}\approx \displaystyle \frac{6{{GM}}^{2}}{5R}\left(1-\displaystyle \frac{1}{45}{e}^{4}-\displaystyle \frac{62}{2835}{e}^{6}\,-\,\cdots \right),\end{eqnarray} \tag{ 36 }$$

such that E_G of the the prolate and oblate is always less than that of the sphere at infinite separation of the two objects. In the same way, it is possible to also show that, for cases (b) and (c), E_G is less than that of the sphere for any value of b.

Zoom In Zoom Out Reset image size

However, as stated above, this is not true for the other cases—it is possible for the value of E_G for the spheroidal configurations (a) and (d) to be greater than that of the sphere. This is further illustrated in figures 5(a) and (b), which are contour plots of ${E}_{G}^{a)}/{E}_{G}^{{\rm{sphere}}}$ and ${E}_{G}^{b)}/{E}_{G}^{{\rm{sphere}}}$ for values of ranging from 0 to 1 (i.e. any ellipticity) and for the displacement b ranging from 0 to 12R. In particular, for small displacements it is possible for the spheroidal E_G to be a factor greater than the spherical case. Taking the oblate spheroid and sphere to have the same volume and mass, then in the limit of $\epsilon \ll 1$ and $b\ll R$, we find, using (23) and (29), that ${E}_{G}^{a)}=3\times {E}_{G}^{{\rm{sphere}}}$ ¹⁴ . Such a factor would already be approximately satisfied when  ≈ 0.01 and d ≈ 0.01R, which could be possible in near-future experiments. For example, $b\ll R$ is satisfied in the proposed nano-sphere experiment [58], such that, taking an oblate spheroid with  ≈ 0.01 rather than a sphere would increase E_G accordingly.

Zoom In Zoom Out Reset image size

The findings here suggest that, for tests of the GQSR process, it may be preferable to use spheroidal rather than spherical mass distributions in nano/micro-object experiments. Figure 5(c) also provides a comparison of ${E}_{G}^{a)}$ against ${E}_{G}^{d)}$, illustrating that, although ${E}_{G}^{a)}$ is generically larger than ${E}_{G}^{d)}$, there is a certain region of parameter space where the opposite is the case. Since a prolate and oblate spheroid can be used to approximate, respectively, a rod or disc for high ellipticity, the self-energy of the difference of these objects could be used to approximate that which could be observed in nano/micro-rod and -disc experiments. Note also that, if an experiment were able to observe state reduction in disagreement with standard QT, then comparing the results for different spheroidal geometries could be used to distinguish the GQSR considered here from other collapse models since we have a direct prediction for how E_G changes with just the ellipticity of the object (e = 0 for a sphere).

In addition to nano/micro-object experiments, it may also be possible to test the GQSR process considered here using BECs. Advantages of these systems include the fact that they are highly controllable systems and have large coherence times due to their extremely low temperatures and high isolation from their environments. Certain superposition states have also already been observed for these experiments, such as a coherent state separated by over 0.5 m, and there are several suggested techniques for creating macroscopic superposition states (see section 3.2).

In section 3.1, we calculate the self-energy of the difference between spherical and spheroidal BEC mass distributions, which are created using harmonic trapping potentials. We then compare, in section 3.3, the corresponding rate of state reduction to the decoherence rate of prominent channels of environmental decoherence in BEC experiments, providing estimates for the values of experimental parameters, such as temperature and scattering length, that would be required to test the GQSR process.

3.1. Gravitational self-energy of the difference between BEC mass distributions

In sections 2.1 and 2.2, we calculated E_G for uniform spherical and spheroidal mass distributions. Although such distributions can be created in nano/micro-object experiments, spherical and spheroidal BEC distributions are generically non-uniform. This non-uniformity is due to the trapping potential that constrains the BEC. For example, to create a spherical BEC, the potential is $V(r)=\tfrac{1}{2}m{\omega }_{0}^{2}{r}^{2}$ where m is the atomic mass, ω₀ is the trapping frequency, and $r=\sqrt{{x}^{2}+{y}^{2}+{z}^{2}}$ is the radial distance from the center of the trap.

Taking the BEC to obey the time-independent Gross–Pitaevskii equation [59, 60]:

$$\begin{eqnarray}&&\left[-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}^{2}+V({\boldsymbol{r}})+{gn}({\boldsymbol{r}})\right]{\psi }_{0}({\boldsymbol{r}})=\mu {\psi }_{0}({\boldsymbol{r}}),\end{eqnarray} \tag{ 37 }$$

where ${\psi }_{0}({\boldsymbol{r}})$ is the BEC wavefunction, μ is the chemical potential of the condensate, g = 4πℏ²a_s/m is the s-wave interaction coupling constant with a_s the s-wave scattering length, and $n({\boldsymbol{r}})=| {\psi }_{0}({\boldsymbol{r}}){| }^{2}$ is the condensate number density; we can solve for the density of the BEC at zero temperature. Here we consider two analytical limiting cases: (1) the Gaussian approximation where we assume that the wavefunction ψ₀ is Gaussian, which is exact for an ideal Bose gas when we neglect the interaction term, and can also be used in describing repulsive BECs with low effective interaction strength, as well as attractive (a_s < 0) BECs [61–66]; and (2) the Thomas–Fermi (TF) approximation [67, 68], which is most appropriate for repulsive BECs (a_s > 0) with large numbers of atoms, where we neglect the kinetic term of (37) in comparison to the interaction term.

3.1.1. Density in the Gaussian approximation

In this section, we consider the Gaussian approximation for BECs, which is useful for characterizing BECs with very low interaction strengths, as well as attractive interactions, as described above. When the interaction term is entirely neglected in the Gross–Pitaevskii equation, we have an ideal Bose gas, and the solution for a general harmonic trapping potential $V({\boldsymbol{r}})=\tfrac{1}{2}m({\omega }_{x}^{2}{x}^{2}+{\omega }_{y}^{2}{y}^{2}+{\omega }_{z}^{2}{z}^{2})$ is:

$$\begin{eqnarray}&&{\psi }_{0}({\boldsymbol{r}})=\sqrt{N}{\left(\displaystyle \frac{m{\omega }_{0}}{{\hslash }\pi }\right)}^{3/4}{{\rm{e}}}^{-\tfrac{1}{2}m\left({\omega }_{x}{x}^{2}+{\omega }_{y}{y}^{2}+{\omega }_{z}{z}^{2}\right)/{\hslash }},\end{eqnarray} \tag{ 38 }$$

where N is the number of condensate atoms, ${\omega }_{0}:= {\left({\omega }_{x}{\omega }_{y}{\omega }_{z}\right)}^{1/3}$ is the geometric average of the trapping frequencies, and the chemical potential is $\mu =\tfrac{1}{2}{\hslash }\left({\omega }_{x}+{\omega }_{y}+{\omega }_{z}\right)$. Taking a spherical trap (ω₀ := ω_x = ω_y = ω_z), the mass density ${\rho }_{0}({\boldsymbol{r}}):= {mn}({\boldsymbol{r}})$ of the condensate is then

$$\begin{eqnarray}&&{\rho }_{0}^{{\rm{sphere}}}({\boldsymbol{r}})=\displaystyle \frac{4}{3\sqrt{\pi }}{\rho }_{0}^{{\rm{sphere}}}{{\rm{e}}}^{-{r}^{2}/{R}_{0}^{2}},\end{eqnarray} \tag{ 39 }$$

where ${\rho }_{0}^{{\rm{sphere}}}:= M/((4/3)\pi {R}_{0}^{3})$, with M = mN the total mass, and

$$\begin{eqnarray}&&{R}_{0}:= \sqrt{{\hslash }/(m{\omega }_{0})}\end{eqnarray} \tag{ 40 }$$

is the width of the Gaussian wavefunction (38). To describe a BEC with attractive forces we can use a variational approach with the ansatz that the ground state is of Gaussian form but we now replace R₀ in (38) with

$$\begin{eqnarray}&&{R}_{0}^{{\prime} }:= {\alpha }_{R}{R}_{0},\end{eqnarray} \tag{ 41 }$$

where α_R is a dimensionless variational parameter which fixes the width of the condensate for a given interaction strength [61–66]. The density for such a spherical BEC will then still be approximated by (39) but with R₀ replaced with ${R}_{0}^{{\prime} }$, where ${R}_{0}^{{\prime} }\lt {R}_{0}$ for a BEC with attractive forces.

To generate a spheroidal BEC, the harmonic trapping potential must be of the form

$$\begin{eqnarray}&&V({\boldsymbol{r}})=\displaystyle \frac{1}{2}m\left[{\omega }_{r}^{2}{r}_{\rho }^{2}+{\omega }_{z}^{2}{z}^{2}\right]\equiv \displaystyle \frac{1}{2}m{\omega }_{r}^{2}\left[{r}_{\rho }^{2}+{\lambda }_{\omega }^{2}{z}^{2}\right],\end{eqnarray} \tag{ 42 }$$

where ${r}_{\rho }=\sqrt{{x}^{2}+{y}^{2}}$ is the radial cylindrical coordinate, ω_r is the radial trapping frequency, ω_z is the axial trapping frequency, and λ_ω := ω_z/ω_r (which is sometimes referred to as the 'asymmetry parameter'). For a prolate spheroid (often called a 'cigar' BEC), λ_ω < 0, whereas, for an oblate spheroid (often called a 'pancake' BEC), λ_ω > 0. In the Gaussian approximation, the density of the BEC is given by

$$\begin{eqnarray}&&{\rho }_{0}^{{\rm{spheroid}}}({\boldsymbol{r}})=\displaystyle \frac{4}{3\sqrt{\pi }}{\rho }_{0}^{{\rm{spheroid}}}{{\rm{e}}}^{-{r}_{\rho }^{2}/{a}_{0}^{{\prime} 2}-{z}^{2}/{c}_{0}^{{\prime} 2}},\end{eqnarray} \tag{ 43 }$$

where ${\rho }_{0}^{{\rm{spheroid}}}:= M/((4/3)\pi {a}_{0}^{{\prime} 2}c{{\prime} }_{0})$, $c{{\prime} }_{0}:= {\alpha }_{c}{c}_{0}$ and $a{{\prime} }_{0}={\alpha }_{a}{a}_{0}$, with ${c}_{0}:= \sqrt{{\hslash }/(m{\omega }_{z})}$ and ${a}_{0}:= \sqrt{{\hslash }/(m{\omega }_{r})}$. Similar to the spherical case, the factors α_a and α_c control the size of the condensate for the given interaction strength.

For a BEC with attractive interactions, the condensate becomes unstable if the number of atoms exceeds a critical value. For a harmonic trap at zero temperature, this critical value can be estimated as [69]

$$\begin{eqnarray}&&{N}_{c}\approx {k}_{c}\displaystyle \frac{{s}_{0}}{| {a}_{s}| },\end{eqnarray} \tag{ 44 }$$

where s₀ is the width of the ground state Gaussian wavefunction of an atom in a parabolic potential well $\sqrt{{\hslash }/(m{\omega }_{0})}$, with ${\omega }_{0}:= {\left({\omega }_{x}{\omega }_{y}{\omega }_{z}\right)}^{1/3}$, and k_c is a constant, which is estimated as ≈0.6 for a single-well spherical trap. See [70] for values of k_c for a double-well trap and, for example [71], for an analytical expression for k_c. Note that (44) is not applicable with the TF approximation since the kinetic term is required to stabilize the system here.

3.1.2. Density in the TF approximation

In this section, we consider the TF approximation for BECs, which is useful for characterizing BECs with strong repulsive interactions, as described in section 3.1 above. Assuming a spherical trap and that the kinetic term can be neglected in comparison to the interaction term, the solution ψ₀ of the Gross–Pitaevskii equation can be used to find the mass density of the BEC in this TF approximation:

$$\begin{eqnarray}&&{\rho }_{0}^{{\rm{sphere}}}({\boldsymbol{r}})=\displaystyle \frac{5}{2}{\rho }_{0}^{{\rm{sphere}}}(1-{r}^{2}/{R}^{2}),\end{eqnarray} \tag{ 45 }$$

where R the radius of the spherical BEC (we assume $\rho ({\boldsymbol{r}})$ is vanishing here), ${\rho }_{0}^{{\rm{sphere}}}:= M/((4/3)\pi {R}^{3})$, and M the total mass of the condensate:

$$\begin{eqnarray}&&M:= {mN}={\int }_{0}^{2\pi }{\int }_{0}^{\pi }{\int }_{0}^{R}\rho ({\boldsymbol{r}}){r}^{2}\sin \theta {\rm{d}}{r}{\rm{d}}\theta {\rm{d}}\phi .\end{eqnarray} \tag{ 46 }$$

In terms of experimental parameters, R is given by

$$\begin{eqnarray}&&R={\left(15{{Na}}_{s}{R}_{0}^{4}\right)}^{1/5},\end{eqnarray} \tag{ 47 }$$

where R₀ is defined in the previous section. The density of the BEC sphere is illustrated in figure 6(a) where contours represent surfaces of constant density, which are spherical surfaces. In contrast to the Gaussian approximation, the Bose gas will clearly be of a larger size in this strongly, repulsive interaction regime.

Zoom In Zoom Out Reset image size

The density function for a spheroidal BEC in the TF approximation can be found by inserting the potential (42) into the Gross–Pitaevskii equation (37), and dropping the kinetic term, to obtain:

$$\begin{eqnarray}&&{\rho }_{0}^{{\rm{spheroid}}}({\boldsymbol{r}})=\displaystyle \frac{5}{2}{\rho }_{0}^{{\rm{spheroid}}}(1-{r}_{\rho }^{2}/{a}^{2}-{z}^{2}/{c}^{2}),\end{eqnarray} \tag{ 48 }$$

where ${\rho }_{0}^{{\rm{spheroid}}}:= M/((4/3)\pi {a}^{2}c)$. This density distribution is illustrated in figures 6(b) and (c) for a prolate (a < c) and oblate (a > c) spheroid where contours are surfaces of constant density, which are similar-shaped spheroidal surfaces. In terms of experimental parameters, the equatorial and polar radii are:

$$\begin{eqnarray}&&a={\left(15{{Na}}_{s}{a}_{0}^{4}{\lambda }_{\omega }\right)}^{1/5},\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}&&c=a/{\lambda }_{\omega },\end{eqnarray} \tag{ 50 }$$

where a₀ is defined in the previous section for the Gaussian approximation. For the TF approximation to be very good, we require that [66]:

$$\begin{eqnarray}&&N\gg \displaystyle \frac{{s}_{0}}{{a}_{s}}.\end{eqnarray} \tag{ 51 }$$

Furthermore, the TF approximation is less accurate near the boundaries of the condensate. Here the density abruptly vanishes in the TF approximation, but in reality there is a more gradual decrease such that the condensate wavefunction will tend to, but never actually reach, zero [72].

3.1.3. Self-energy of the difference between spherical BECs

Now that we have mass distributions for BEC spheres and spheroids, we can determine the value of E_G for the different shapes and density functions using (18) or (19). An approach to this is discussed in appendices B–F where we also calculate the gravitational potential of these objects. For the spherical BEC in the Gaussian approximation (see appendix C), we find

$$\begin{eqnarray}&&{E}_{G}=\displaystyle \frac{{{GM}}^{2}}{{R}_{0}^{{\prime} }}\left(\sqrt{\displaystyle \frac{2}{\pi }}-\displaystyle \frac{1}{2{\lambda }_{0}}\mathrm{erf}(\sqrt{2}{\lambda }_{0})\right),\end{eqnarray} \tag{ 52 }$$

where we have defined λ₀ := b/(2R'₀) with ${R}_{0}^{{\prime} }$ given by (41). In contrast, in the TF regime, we obtain (see appendix B):

$$\begin{eqnarray}{E}_{G}=\left\{\begin{array}{cl}\displaystyle \frac{10{{GM}}^{2}}{7R}\left(2{\lambda }^{2}-\displaystyle \frac{21}{5}{\lambda }^{4}+\displaystyle \frac{7}{2}{\lambda }^{5}-\displaystyle \frac{3}{4}{\lambda }^{7}+\displaystyle \frac{1}{10}{\lambda }^{9}\right) & \mathrm{if}\ 0\leqslant \lambda \leqslant 1,\\ \displaystyle \frac{10{{GM}}^{2}}{7R}\left(1-\displaystyle \frac{7}{20\lambda }\right) & \mathrm{if}\ \lambda \geqslant 1,\end{array}\right.\end{eqnarray} \tag{ 53 }$$

where λ := b/(2R) and R is given by (47).

These self-energy differences are illustrated in figure 7. For the same total mass and volume (and so average density), E_G of a spherical BEC in the TF regime is always greater than that of a uniform one. This is exemplified by figure 7(a) and is due to the fact that the density is more constrained towards the center. The fact that E_G is different despite the potential outside a non-uniform sphere being the same as a uniform sphere, could provide a further possibility for distinguishing the state reduction process considered here to other models.

Zoom In Zoom Out Reset image size

In figure 7(b), we plot E_G of a spherical BEC in the TF regime against the Gaussian regime for a ¹³³Cs BEC with 10⁶ atoms, the same trapping frequency ${\omega }_{0}\,=\,100\,\mathrm{Hz}$, and with the standard scattering length in the former regime, but with zero scattering in the latter so that we have an ideal BEC. For all values of b, the Gaussian E_G is always greater than the TF case. This is principally due to R₀ being much smaller than R in this case, with the gap increasing as N increases. Therefore, with attractive interactions, we would expect E_G to increase further under the condition that all other BEC parameters, apart from the scattering length, stay the same.

3.1.4. Self-energy of the difference between spheroidal BECs

The generic value of E_G for spheroidal BECs is more complicated (see appendices E and F) and here we just provide the expression for configuration (b) in figure 3 (the symmetry-axis displaced prolate) for the TF regime in the limit of high ellipticity (i.e. to first order in , which is defined by (24)):

$$\begin{eqnarray}{E}_{G}=\left\{\begin{array}{cl}\displaystyle \frac{10{{GM}}^{2}}{7a}\left({\mathtt{A}}/24-{\mathtt{B}}\mathrm{ln}\epsilon \right) & \mathrm{if}\ 0\leqslant \lambda \leqslant 1,\\ \displaystyle \frac{10{{GM}}^{2}}{7a}\left({\mathtt{C}}/1536-\mathrm{ln}\epsilon \right) & \mathrm{if}\ \lambda \geqslant 1,\end{array}\right.\end{eqnarray} \tag{ 54 }$$

where λ := b/(2c) and

$$\begin{eqnarray}&&{\mathtt{A}}:= 144{\lambda }^{2}(\mathrm{ln}2-1)-168{\lambda }^{4}(3\mathrm{ln}2-4)-378{\lambda }^{5}+133{\lambda }^{7}-25{\lambda }^{9},\end{eqnarray} \tag{ 55 }$$

$$\begin{eqnarray}&&{\mathtt{B}}:= 10\left(6{\lambda }^{2}-21{\lambda }^{4}+21{\lambda }^{6}-6{\lambda }^{7}+{\lambda }^{9}\right),\end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\mathtt{C}} & := & 6\gamma {\coth }^{-1}(1-2\lambda )-6\kappa {\coth }^{-1}(1+2\lambda )\\ & & -768\left(1-6{\lambda }^{2}+21{\lambda }^{4}\right){\coth }^{-1}\lambda \\ & & -64\left(12{\lambda }^{7}-66{\lambda }^{5}-284{\lambda }^{3}+81\lambda -24\mathrm{ln}2\right)\\ & & -315\lambda (16{\lambda }^{2}-3)\mathrm{ln}\displaystyle \frac{{\lambda }^{2}-1}{{\lambda }^{2}}\end{array}\end{eqnarray} \tag{ 57 }$$

with

$$\begin{eqnarray}&&\gamma := {\left(1-2\lambda \right)}^{3}\left[128+\lambda \left(2\lambda \left[207+4\lambda \left(2\lambda (\lambda +3)(2\lambda -3)\right)\right]+453\right)\right],\end{eqnarray} \tag{ 58 }$$

$$\begin{eqnarray}&&\kappa := {\left(1+2\lambda \right)}^{3}\left[128+\lambda \left(2\lambda [207+4\lambda (2\lambda (\lambda -3)(2\lambda +3))]-453\right)\right].\end{eqnarray} \tag{ 59 }$$

Assuming the TF regime, the value of E_G for the four configurations (a)–(d) (see figure 3) of BEC spheroids is plotted in figures 8(a) and (b) against the BEC sphere case for  = 0.5 and  = 0.01. As in the uniform case, the value of E_G for configurations (b) and (c) is always less than that of a BEC sphere, whereas, the other spheroidal configurations can have larger E_G values at certain displacements and ellipticity values. However, again, the sphere always has the greatest E_G at infinity—the values of E_G in the BEC TF case compared to the uniform case (32)–(34) are just 25/21 ≈ 1.2 larger for each object.

Zoom In Zoom Out Reset image size

In figures 9(a)–(c) we compare E_G of the BEC sphere with a spheroidal BEC in configuration (a), the BEC sphere with a spheroidal BEC in configuration (d), and the spheroidal configuration (a) with (d), for all BECs in the TF regime and assuming the same volume and density for the different objects. These are very similar to the uniform cases 5(a), (b) and illustrate again that it may be preferable to use spheroidal rather than spherical objects for testing GQSR.

Zoom In Zoom Out Reset image size

In figure 10, we also plot spherical and spheroidal configurations (a) and (b) for BECs in the Gaussian approximation with  = 0.75 (e ≈ 0.7) and displacement b from zero to 10R. As in the TF approximation, the oblate case can have a value of E_G that is greater than the spherical case. Note that for high values of ellipticity, it is possible to enter a quasi-one or two dimensional regime where the quantum and thermal motion can be frozen in two or one dimensions (see e.g. [73]). This is to be distinguished from the case when the BEC looks lower dimensional from only a geometrical point of view but locally has a three-dimensional character. In certain configurations, it can be a good approximation to neglect the spatial dependence of the density in one or two dimensions, potentially simplifying the calculation of E_G for such BEC states.

Zoom In Zoom Out Reset image size

3.1.5. Self-energy difference in BEC experiments

Now that we have calculated E_G for mass distributions that can be generated by BEC experiments, let us consider what sort of experimental parameters might be required to test the gravitationally-induced state reduction model. Taking a spherical BEC for simplicity, when the separation of the two BEC states is of the order of their (effective) diameter, the value of E_G is of the order (assuming γ = 1/(8π) in (18)):

$$\begin{eqnarray}&&{E}_{G}\sim \displaystyle \frac{{{Gm}}^{2}{N}^{2}}{R}.\end{eqnarray} \tag{ 60 }$$

For example, in the TF approximation, when two spherical BEC states are touching, the value of E_G is found to be:

$$\begin{eqnarray}&&{E}_{G}=\displaystyle \frac{13{{Gm}}^{2}{N}^{2}}{14R}.\end{eqnarray} \tag{ 61 }$$

Using this expression for E_G, for a ¹³³Cs BEC of radius 1 μm, we would need around 4 × 10⁹ atoms in each state for a collapse lifetime of around 2 s. In (60), there is a stronger dependence on the number of atoms N than on the radius R, and so N = 4 × 10¹⁰ and $R=0.1\,\mathrm{mm}$ would also cause the same collapse rate, while potentially being more experimentally feasible due to the reduced density (see section 3.3). On the other hand, if, for example, γ = 8π were found to be more appropriate in (18), then a collapse time of $2\,{\rm{s}}$ would occur when N ≈ 10⁹ and $R=0.1\,\mathrm{mm}$ or N ≈ 10⁸ and $R=1\,\mu {\rm{m}}$. Allowing for smaller timescales than $2\,{\rm{s}}$ would also improve the required values for N and R.

Although such numbers of atoms have not been achieved yet for a ¹³³Cs BEC experiment, over 10⁹ atoms were reported for a hydrogen BEC in 1998 [74], and over 10⁸ atoms for a ²³Na BEC in 2006 [75] (also see [76] for a ²³Na BEC of over 10⁷ atoms in 2004). These were single-well rather than double-well BECs, and so not large macroscopic superposition states. However, in sections 3.2 and 3.3 we discuss how large macroscopic states, such as NOON states, or approximations to these, could be generated in double-well BECs, and what sort of experimental parameters would be required in order for GQSR to be seen in the presence of environmental decoherence.

3.2. Generating macroscopic superposition states with double-well BECs

A double-well BEC can, in principle, be used to create a macroscopic superposition state (see, for example, [37–40, 77–83]). The full Hamiltonian of this system is:

$$\begin{eqnarray}&&\hat{H}=\int {{\rm{d}}}^{3}{\boldsymbol{r}}{\hat{{\rm{\Psi }}}}^{\dagger }({\boldsymbol{r}}){\hat{H}}_{{\rm{DW}}}\hat{{\rm{\Psi }}}({\boldsymbol{r}})+\displaystyle \frac{1}{2}g\int {{\rm{d}}}^{3}{\boldsymbol{r}}{\hat{{\rm{\Psi }}}}^{\dagger }({\boldsymbol{r}}){\hat{{\rm{\Psi }}}}^{\dagger }({\boldsymbol{r}})\hat{{\rm{\Psi }}}({\boldsymbol{r}})\hat{{\rm{\Psi }}}({\boldsymbol{r}}),\end{eqnarray} \tag{ 62 }$$

where

$$\begin{eqnarray}&&\left[\hat{{\rm{\Psi }}}({\boldsymbol{r}}),{\hat{{\rm{\Psi }}}}^{\dagger }({\boldsymbol{r}}^{\prime} )\right]=\delta ({\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} );\end{eqnarray} \tag{ 63 }$$

$$\begin{eqnarray}&&{\hat{H}}_{{\rm{DW}}}:= \left[-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}^{2}+{V}_{{\rm{DW}}}({\boldsymbol{r}})\right];\end{eqnarray} \tag{ 64 }$$

V_DW is the particular double-well potential, which we take to be symmetric; and we have assumed that the inter-atomic interaction can be well-approximated by two-body s-wave scattering.

Assuming that the energy barrier between the two wells is large enough, we make the ansatz that the BEC can be described as consisting of atoms that occupy a condensed state $| {\psi }_{L}\rangle$ of the left well, or a condensed state $| {\psi }_{R}\rangle$ of the right well, which are taken to be approximately orthogonal, $\langle {\psi }_{L}| {\psi }_{R}\rangle \approx 0$. That is, we assume that $\hat{{\rm{\Psi }}}$ can be approximated by:

$$\begin{eqnarray}&&\hat{{\rm{\Psi }}}({\boldsymbol{r}},t)={\psi }_{L}({\boldsymbol{r}},t){\hat{a}}_{L}(t)+{\psi }_{R}({\boldsymbol{r}},t){\hat{a}}_{R}(t),\end{eqnarray} \tag{ 65 }$$

where ${\hat{a}}_{L}$ and ${\hat{a}}_{R}$ are the annihilation operators for the states $| \psi {\rangle }_{L}$ and $| \psi {\rangle }_{R}$, which have localized wavefunctions ψ_L and ψ_R. These obey the usual bosonic commutation rules:

$$\begin{eqnarray}&&[{\hat{a}}_{i},{\hat{a}}_{j}^{\dagger }]={\delta }_{{ij}},\hspace{2cm}[{\hat{a}}_{i},{\hat{a}}_{j}]=0,\end{eqnarray} \tag{ 66 }$$

where $i,j\,=\,L,R$. The wavefunctions ψ_L and ψ_R are assumed to have negligible overlap such that they are approximately orthogonal:

$$\begin{eqnarray}&&\int {\rm{d}}{\boldsymbol{r}}{\psi }_{i}^{* }({\boldsymbol{r}},t){\psi }_{j}({\boldsymbol{r}},t)\approx {\delta }_{{ij}},\end{eqnarray} \tag{ 67 }$$

and the operator for the total number of particles, which is conserved, is given by

$$\begin{eqnarray}&&\hat{N}:= \int {\rm{d}}{\boldsymbol{r}}\hat{{\rm{\Psi }}}{\hat{{\rm{\Psi }}}}^{\dagger }={\hat{a}}_{L}^{\dagger }{\hat{a}}_{L}+{\hat{a}}_{R}^{\dagger }{\hat{a}}_{R}.\end{eqnarray} \tag{ 68 }$$

In the nonlinear tight-binding approximation [84], an adiabatic approximation is applied where ψ_L and ψ_R are real and their spatial profiles adapt adiabatically to the instantaneous number of particles. In this tight-binding approximation, the wavefunctions depend implicitly on time t through the number of particles in each well ${N}_{i}=\langle {\hat{a}}_{i}^{\dagger }{\hat{a}}_{i}\rangle$, i.e. ${\psi }_{i}({\boldsymbol{r}},t)={\psi }_{i}({\boldsymbol{r}},{N}_{i}(t))$. In our large separation approximation, and assuming macroscopic occupation of the two states, the wavefunctions ψ_L and ψ_R (when multiplied by $\sqrt{{N}_{L}}$ and $\sqrt{{N}_{R}}$) obey, to a good approximation, solutions of the Gross–Pitaevskii equation (37) with potential V_DW [84] (for, alternatively, a full variational approach, see [85])¹⁵ .

Plugging our ansatz (65) into the Hamiltonian (62), we obtain

$$\begin{eqnarray}&&\hat{H}={\hat{H}}_{1}+{\hat{H}}_{2},\end{eqnarray} \tag{ 69 }$$

where

$$\begin{eqnarray}&&{\hat{H}}_{1}:= {\xi }_{L}{\hat{a}}_{L}^{\dagger }{\hat{a}}_{L}+{\xi }_{R}{\hat{a}}_{R}^{\dagger }{\hat{a}}_{R}+{J}_{{LR}}({\hat{a}}_{L}^{\dagger }{\hat{a}}_{R}+{\hat{a}}_{R}^{\dagger }{\hat{a}}_{L}),\end{eqnarray} \tag{ 70 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\hat{H}}_{2} & := & {U}_{L}{\hat{a}}_{L}^{\dagger 2}{\hat{a}}_{L}^{2}+{U}_{R}{\hat{a}}_{R}^{\dagger 2}{\hat{a}}_{R}^{2}+4{U}_{{LRLR}}{\hat{a}}_{L}^{\dagger }{\hat{a}}_{R}^{\dagger }{\hat{a}}_{L}{\hat{a}}_{R}\\ & & +2{U}_{{LLLR}}\left({\hat{a}}_{L}^{\dagger 2}{\hat{a}}_{L}{\hat{a}}_{R}+\,{\rm{h}}.{\rm{c}}.\right)+2{U}_{{RRRL}}\left({\hat{a}}_{R}^{\dagger 2}{\hat{a}}_{R}{\hat{a}}_{L}+\,{\rm{h}}.{\rm{c}}.\right)\\ & & +{U}_{{LLRR}}\left({\hat{a}}_{L}^{\dagger 2}{\hat{a}}_{R}^{\dagger 2}+\,{\rm{h}}.{\rm{c}}.\right),\end{array}\end{eqnarray} \tag{ 71 }$$

with

$$\begin{eqnarray}&&{\xi }_{L}:= \int {{\rm{d}}}^{3}{\boldsymbol{r}}{\psi }_{L}{\hat{H}}_{{\rm{DW}}}{\psi }_{L},{\xi }_{R}:= \int {{\rm{d}}}^{3}{\boldsymbol{r}}{\psi }_{R}{\hat{H}}_{{\rm{DW}}}{\psi }_{R},\end{eqnarray} \tag{ 72 }$$

$$\begin{eqnarray}&&{J}_{{LR}}:= \int {{\rm{d}}}^{3}{\boldsymbol{r}}{\psi }_{L}{\hat{H}}_{{\rm{DW}}}{\psi }_{R},{U}_{L}:= g\int {{\rm{d}}}^{3}{\boldsymbol{r}}{\psi }_{L}^{4},{U}_{R}:= g\int {{\rm{d}}}^{3}{\boldsymbol{r}}{\psi }_{R}^{4},\end{eqnarray} \tag{ 73 }$$

$$\begin{eqnarray}&&{U}_{{LRLR}}:= g\int {{\rm{d}}}^{3}{\boldsymbol{r}}{\psi }_{L}^{2}{\psi }_{R}^{2},{U}_{{LLRR}}:= g\int {{\rm{d}}}^{3}{\boldsymbol{r}}{\psi }_{L}^{2}{\psi }_{R}^{2},\end{eqnarray} \tag{ 74 }$$

$$\begin{eqnarray}&&{U}_{{LLLR}}:= g\int {{\rm{d}}}^{3}{\boldsymbol{r}}{\psi }_{L}^{2}{\psi }_{L}{\psi }_{R},{U}_{{RRRL}}:= g\int {{\rm{d}}}^{3}{\boldsymbol{r}}{\psi }_{R}^{2}{\psi }_{R}{\psi }_{L},\end{eqnarray} \tag{ 75 }$$

and we are taking ψ_L and ψ_R to be real. The Hamiltonian (69) can be shown to contain an analytic solution [86]. It can also be approximated by an extended two-mode Bose–Hubbard model in the nonlinear tight-binding approximation [84]. In the case that the spatial profile of ψ_L and ψ_R is approximately independent of the number of particles in each well (the standard tight-binding approximation), it can be further approximated with the two-mode version of the Bose–Hubbard model [73, 84, 87–89]:

$$\begin{eqnarray}&&\hat{H}={E}_{{LR}}({\hat{a}}_{L}^{\dagger }{\hat{a}}_{R}+{\hat{a}}_{R}^{\dagger }{\hat{a}}_{L})+\displaystyle \frac{1}{2}U({\hat{a}}_{L}^{\dagger 2}{\hat{a}}_{L}^{2}+{\hat{a}}_{R}^{\dagger 2}{\hat{a}}_{R}^{2}),\end{eqnarray} \tag{ 76 }$$

where U = U_L; we have assumed U_L = U_R; we have removed terms proportional to the number operator $\hat{N}$ since this commutes with $\hat{H};$ and we have neglected any atomic collisions in the overlapping region of the two modes. Here the J_LR terms are responsible for quantum tunneling between the two wells, and the U terms are the atom-atom interactions within each well.

There have been several proposals for generating a macroscopic superposition state (Schrödinger cat state) in a double-well BEC. For example, in [37, 38] it is considered that, starting from a repulsive BEC, if the interaction strength g is varied adiabatically to a negative value (using a Feshbach resonance), then a cat state can be prepared. This occurs because a NOON state is the ground state of the two-mode Bose–Hubbard model (76) with strong attractive interactions. In [39], it is shown that the ground state becomes degenerate with the first excited state in this case, such that there needs to be an exponentially long time to create an exact NOON state. However, in [90] it was found that, for realistic parameters and time-scales, an approximate NOON state can be generated with a smooth change in the scattering length. An alternative to this method is to use a Feshbach resonance to make a sudden change in the scattering length [39]. For example, a repulsive BEC could be prepared in a single-well and then the tunneling barrier is raised adiabatically to divide the well into two equal parts (forming a so-called 'coherent' state [73] when neglecting interactions), then a Feshbach resonance is used to suddenly switch g from a positive to a negative value such that the state dynamically evolves to a large macroscopic superposition state.

Another possibility would be to set the scattering length to zero and drive the system to an upper excited state, then slowly increase the interactions (keeping them repulsive) while, at the same time, decreasing the inter-well tunneling to zero [40]. This method is possible since a NOON state is the upper energy state of the repulsive Bose–Hubbard model, and has the advantage that the BEC does not need to move to an attractive state, which can become unstable [66]. Rather than modifying the scattering length, a cat state could also be generated by manipulating the BEC with an external laser [77, 80, 81]. For example, in [81], it is suggested that a far off-resonance laser could be used to imprint a π-phase on one of the wells such that the quantum wavepacket bifurcates. The tunneling barrier is then raised to halt the evolution and fix the cat state.

Once a macroscopic superposition state, such as a NOON state, has been prepared, we need to make sure that we can experimentally distinguish it from a classical statistical mixture. For a double-well BEC, one possibility is to look for a a non-zero Nth-order correlation $\langle {\hat{a}}_{L}^{\dagger N}{\hat{a}}_{R}^{N}\rangle$ [90, 91]. For an exact NOON state, $| \mathrm{NOON}\rangle$, where we have a superposition of N particles in the left-hand state $| {\psi }_{L}\rangle$ and N particles in the right-hand state $| {\psi }_{R}\rangle$, which we write as $(| {N}0\rangle +| 0{N}\rangle )/\sqrt{2}$, we have

$$\begin{eqnarray}&&\langle \mathrm{NOON}| {\hat{a}}_{L}^{\dagger N}{\hat{a}}_{R}^{N}| \mathrm{NOON}\rangle =\displaystyle \frac{N!}{2},\end{eqnarray} \tag{ 77 }$$

whereas, for a statistical mixture, we have zero. Experimental methods for measuring $\langle {\hat{a}}_{L}^{\dagger N}{\hat{a}}_{R}^{N}\rangle$ in double-well BECs can be found in [90, 91].

As well as being able to distinguish a NOON state from a statistical mixture, we also need need to make sure that we can experimentally distinguish the GQSR process from environmental decoherence. That is, we would ideally like the objective collapse rate to be greater than the rate of environmental decoherence. Given an initial NOON state, $| \mathrm{NOON}\rangle$, we can use (21), (22) to write down the density operator for the state under the process of GQSR:

$$\begin{eqnarray}&&\hat{\rho }(t)=\displaystyle \frac{1}{2}{{\rm{e}}}^{-{E}_{G}t/{\hslash }}\left[| {N}0\rangle +| 0{N}\rangle \right]\left[\langle {N}0| +\langle 0{N}| \right]\end{eqnarray} \tag{ 78 }$$

$$\begin{eqnarray}&&+\,\displaystyle \frac{1}{2}\left(1-{{\rm{e}}}^{-{E}_{G}t/{\hslash }}\right)\left[| {N}0\rangle \langle {N}0| +| 0{N}\rangle \langle 0{N}| \right].\end{eqnarray} \tag{ 79 }$$

In terms of the annihilation operators of the left and right states, ${\hat{a}}_{L}$ and ${\hat{a}}_{R}$, we have

$$\begin{eqnarray}&&| \mathrm{NOON}\rangle =\displaystyle \frac{1}{\sqrt{2N!}}({\hat{a}}_{L}^{\dagger N}+{\hat{a}}_{R}^{\dagger N})| 0\rangle ,\end{eqnarray} \tag{ 80 }$$

such that the N-particle correlation $\langle {\hat{a}}_{L}^{\dagger N}{\hat{a}}_{R}^{N}\rangle$ evolves in time as

$$\begin{eqnarray}&&\langle {\hat{a}}_{L}^{\dagger N}{\hat{a}}_{R}^{N}\rangle (t)={{\rm{e}}}^{-{E}_{G}t/{\hslash }}\langle {\hat{a}}_{L}^{\dagger N}{\hat{a}}_{R}^{N}{\rangle }_{0},\end{eqnarray} \tag{ 81 }$$

where $\langle {\hat{a}}_{L}^{\dagger N}{\hat{a}}_{R}^{N}{\rangle }_{0}$ is given by (77). We now compare this evolution of the N-particle correlation to that imposed by various environmental decoherence channels in double-well BECs.

3.3. Environmental decoherence

There are several channels of environmental decoherence in BEC systems. Here we concentrate on the prominent ones due to three-body recombination, interactions with the thermal cloud, and interactions with foreign atoms. We also briefly discuss noise due to the trapping potential.

3.3.1. Three-body recombination

Three-body recombination is the process where three atoms in the condensate collide to form a molecule (atom-atom bound state) and a single atom, which can both then escape the trap. This process often limits the lifetime and size of condensates. In [92], a master equation was derived for three-body loss in the Born–Markov approximation and for a BEC with repulsive interactions. Since this is a three-body problem, this master equation is of the following form for a double-well BEC in the two-mode approximation [93]:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\hat{\rho }(t)}{{\rm{d}}{t}}=-\displaystyle \frac{{\rm{i}}}{{\hslash }}[\hat{H},\hat{\rho }(t)]+{\gamma }_{3}\displaystyle \sum _{k=L,R}[{\hat{a}}_{k}^{3}\hat{\rho }(t){\hat{a}}_{k}^{\dagger 3}-\displaystyle \frac{1}{2}\{{\hat{a}}_{k}^{\dagger 3}{\hat{a}}_{k}^{3},\hat{\rho }(t)\}],\end{eqnarray} \tag{ 82 }$$

where

$$\begin{eqnarray}&&{\gamma }_{3}:= \displaystyle \frac{{K}_{3}}{72}\int {{\rm{d}}}^{3}{\boldsymbol{x}}| {\psi }_{L,R}({\boldsymbol{r}}){| }^{6}\end{eqnarray} \tag{ 83 }$$

$$\begin{eqnarray}&&=\displaystyle \frac{{K}_{3}}{72}{n}^{2},\end{eqnarray} \tag{ 84 }$$

with n the condensate number density and K₃ the recombination event rate, which can be approximated as [94]:

$$\begin{eqnarray}&&{K}_{3}=23\displaystyle \frac{{\hslash }}{m}{a}_{s}^{4}.\end{eqnarray} \tag{ 85 }$$

The N-particle correlation for a NOON state under this master equation is then [93]:

$$\begin{eqnarray}&&\langle {\hat{a}}_{L}^{\dagger N}{\hat{a}}_{R}^{N}\rangle (t)={{\rm{e}}}^{-{\gamma }_{3}{Nt}}\langle {\hat{a}}_{L}^{\dagger N}{\hat{a}}_{R}^{N}{\rangle }_{0}.\end{eqnarray} \tag{ 86 }$$

Comparing to the gravitationally-induced collapse rate for a NOON state (81), we require that

$$\begin{eqnarray}&&{E}_{G}/{\hslash }\gg {\gamma }_{3}N.\end{eqnarray} \tag{ 87 }$$

Taking E_G to be of the form (60) for simplicity, we need

$$\begin{eqnarray}&&\displaystyle \frac{{{Gm}}^{2}N}{R}\gg \displaystyle \frac{23}{72m}{a}_{s}^{4}{{\hslash }}^{2}{n}^{2}.\end{eqnarray} \tag{ 88 }$$

Assuming, for example, a ¹³³Cs BEC with N ∼ 4 × 10⁹ and R ∼ 10 μm (such that τ ∼ 10 s), then to obtain a three-body recombination rate that is ten times slower in the TF regime, we would need to utilize a Feshbach resonance in order to reduce the scattering length by approximately four orders of magnitude (and we take the trapping frequency to be around $300\,\mathrm{Hz}$). Increasing the number of atoms to 4 × 10¹⁰ instead, then a radius $R\sim 0.1\,\mathrm{mm}$ (trap frequency $10\,\mathrm{Hz}$) and a reduction in the scattering length by three orders of magnitude would be enough. Assuming instead a Gaussian approximation, then in order to operate in this regime, the trapping frequency and/or scattering length need to be reduced further, which only lowers the decoherence rate. Note that, if it were found to be more appropriate to take γ = 8π rather than γ = 1/(8π) in (18), then this would increase E_G by almost three orders of magnitude, significantly improving the experimental feasibility. For example, in this case we could have around 6 × 10⁸ atoms and $R=0.1\,\mathrm{mm}$, with the interaction strength reduced by two orders of magnitude.

In several of the proposals to create a NOON state that were discussed in section 3.2, the (attractive or repulsive) interaction strength is modified. For example, in [40], the repulsive interaction strength is increased while the inter-well tunneling is reduced to zero. In this case, once the NOON state is prepared (or a good approximation to one) the interaction strength would then likely have to be reduced in order to prolong the coherence of, at least an approximation to, the state in light of three-body interactions. Alternatively, other methods could be employed to inhibit three-body decay, such as using an external laser [95, 96] or lowering the effective dimensionality of the BEC [97–99]. Here we have assumed a three dimensional BEC throughout. However, although condensation cannot occur in one or two dimensional uniform systems, with a harmonic trap it is possible to have condensation in an ideal Bose gas in two dimensions, and macroscopic occupation of the lowest energy state in one dimension at finite temperatures [66]. These lower dimensional systems can be achieved when one or two of the harmonic trapping frequencies are much higher than the others, i.e. in the limit of a very flat oblate spheroid or thin prolate spheroid. Unlike in three dimensions, in a Bose gas of one or two dimensions, the three-body decay can become temperature dependent and vanishing at absolute zero. Therefore, reducing the effective dimensionality of the trap, and operating at low temperatures may be another possibility for inhibiting decoherence due to three-body decay.

As stated in section 3.1.1, for a BEC with attractive interactions, the condensate becomes unstable if the number of atoms exceeds a critical value N_c, which for a spherical trap at zero temperature is given by (44). Therefore, if a NOON state is formed with an attractive BEC, the number of atoms N needs to be lower than N_c [73]. One possibility is to increase N_c by lowering the scattering length. However, an exact NOON state from this method is only obtained in the limit of infinite attractive interactions [90]. Therefore, a lower a_s would likely lead to an approximation to a NOON state, for which we would have to calculate the rate of GQSR, and will be the concern of future work. It may instead be preferable to utilize one of the methods outlined in section 3.3 that generates a NOON state with a repulsive BEC.

3.3.2. Thermal cloud interactions

Interactions between the condensate atoms and atoms in the thermal cloud (the noncondensed atoms due to a finite temperature) will also lead to decoherence of a NOON state [79, 100]. These interactions can be of three types: single particle loss $C+{NC}\to {NC}+{NC}$, two particle loss $C+C\to {NC}+{NC}$, and scattering $C+{NC}\to C+{NC}$ (together with the opposite processes) [101]. In [100], assuming a Born–Markov and standard tight-binding approximation, a quantum master equation was derived for the scattering process, where the thermal cloud environment can learn the quantum state of the condensate system. This is of the form:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\hat{\rho }(t)}{{\rm{d}}{t}}=-\displaystyle \frac{{\rm{i}}}{{\hslash }}\left[\hat{H},\hat{\rho }(t)\right]-{\gamma }_{t}\left[{\hat{a}}_{L}^{\dagger }{\hat{a}}_{L}-{\hat{a}}_{R}^{\dagger }{\hat{a}}_{R},\left[{\hat{a}}_{L}^{\dagger }{\hat{a}}_{L}-{\hat{a}}_{R}^{\dagger }{\hat{a}}_{R},\hat{\rho }(t)\right]\right],\end{eqnarray} \tag{ 89 }$$

where

$$\begin{eqnarray}&&{\gamma }_{t}:= 64{\pi }^{4}{a}_{s}^{2}{n}_{{\rm{th}}}{v}_{t},\end{eqnarray} \tag{ 90 }$$

with ${v}_{t}:= \sqrt{2{k}_{B}T/m}$ being the thermal velocity of the atoms in the thermal cloud; T the temperature; and n_th the thermal cloud number density, which can be approximated by:

$$\begin{eqnarray}&&{n}_{{\rm{th}}}=\displaystyle \frac{{{\rm{e}}}^{-\mu /{k}_{B}T}}{{V}_{{\rm{th}}}}{\left(\displaystyle \frac{{k}_{B}T}{{\hslash }\omega }\right)}^{3},\end{eqnarray} \tag{ 91 }$$

$$\begin{eqnarray}&&{V}_{{\rm{th}}}:= \displaystyle \frac{4}{3}\pi {R}_{{\rm{th}}}^{3},\end{eqnarray} \tag{ 92 }$$

$$\begin{eqnarray}&&{R}_{{\rm{th}}}:= \sqrt{\displaystyle \frac{2{k}_{B}T}{m{\omega }^{2}}},\end{eqnarray} \tag{ 93 }$$

where $\omega := {\left({\omega }_{x}{\omega }_{y}{\omega }_{z}\right)}^{1/3}$ and ω_{x, y, z} are the various harmonic trapping frequencies. In [102, 103], a master equation was derived for the scattering of thermal particles off a single condensate within the TF regime. The rate γ_t in this case becomes [103]:

$$\begin{eqnarray}&&{\gamma }_{t}=\displaystyle \frac{4{k}_{B}T{\mu }_{{\rm{TF}}}^{4}}{9{\pi }^{4}{{\hslash }}^{5}{\omega }^{4}{N}^{2}}{{\rm{e}}}^{\mu /{k}_{B}T},\end{eqnarray} \tag{ 94 }$$

where μ is the chemical potential of the non-condensed cloud, we have assume a spherical trapping potential, and

$$\begin{eqnarray}&&{\mu }_{{\rm{TF}}}=\displaystyle \frac{1}{2}{\hslash }\omega {\left(15{Na}/{R}_{0}\right)}^{2/5}\end{eqnarray} \tag{ 95 }$$

is the chemical potential of a spherical BEC in the TF regime.

The N-particle correlation for a NOON state under the master equation (89) is:

$$\begin{eqnarray}&&\langle {\hat{a}}_{L}^{\dagger N}{\hat{a}}_{R}^{N}\rangle (t)={{\rm{e}}}^{-{\gamma }_{t}{N}^{2}t}\langle {\hat{a}}_{L}^{\dagger N}{\hat{a}}_{R}^{N}{\rangle }_{0}.\end{eqnarray} \tag{ 96 }$$

Taking E_G to be of the form (60), we require

$$\begin{eqnarray}&&\displaystyle \frac{{{Gm}}^{2}}{R{\hslash }}\gg {\gamma }_{t}.\end{eqnarray} \tag{ 97 }$$

Assuming a Gaussian ¹³³Cs BEC with the interaction strength reduced by six orders of magnitude and 4 × 10⁹ atoms, as considered as a possibility in the previous section, we would need to increase the trapping potential so that the radius is of order $1\,\mu {\rm{m}}$, and operate at a temperature $T\lesssim 1\,\mathrm{nK}$. A temperature of $0.5\,\mathrm{nK}$ has been achieved for a low-density ²³Na BEC in a single-well potential [104]. If, on the other hand, we want to work in the TF regime, then more challenging experimental parameters appear to be necessary. For example, environmental decoherence would be five times slower than collapse when there is around 4 × 10¹¹ atoms in a condensate of radius of $0.1\,\mathrm{mm}$, the interaction strength is reduced by a further two orders of magnitude as compared to the TF regime considered in the previous section, and the temperature is $T\lesssim 0.1\,\mathrm{nK}$. Therefore, as suggested in the previous section, if a NOON state is prepared by changing the interaction strength then, to prolong the lifetime of the state, it would be preferable to subsequently significantly reduce the interaction strength so that we are working with an approximately ideal gas.

The temperature bound and/or the condensate radius can be increased in the Gaussian approximation by further lowering the interaction strength and, at the same time, either keeping the total atom number the same or increasing it. Also, trap engineering and symmetrization of the environment would help [100]. However, it is possible that a Born–Markov approximation is not appropriate for the description of thermal cloud decoherence in this case, such that the estimates provided here would be inaccurate [39]. Furthermore, as discussed in the previous section, it is possible that these values could be improved if we took γ = 8π rather than γ = 1/(8π) in (18). For example, in this case it would be possible to lower the total atom number to 10⁸ while keeping the rest of the parameters the same.

3.3.3. Foreign atom interactions

Decoherence can also occur due to interactions with background gas particles at room temperature¹⁶ . These foreign particles collide with the condensate atoms and can either cause them to leave the trap entirely or heat up [105]. Assuming that all collisions cause atoms to leave the condensate, a master equation for this process was derived in [106] assuming a Born–Markov approximation. Since this is a one-body process, this master equation is of the form:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\hat{\rho }(t)}{{\rm{d}}{t}}=-\displaystyle \frac{{\rm{i}}}{{\hslash }}[\hat{H},\hat{\rho }(t)]+{\gamma }_{f}\displaystyle \sum _{k=L,R}[{\hat{a}}_{k}\hat{\rho }(t){\hat{a}}_{k}^{\dagger }-\displaystyle \frac{1}{2}\{{\hat{a}}_{k}^{\dagger }{\hat{a}}_{k},\hat{\rho }(t)\}].\end{eqnarray} \tag{ 98 }$$

A rough estimate of the rate γ_f can be calculated assuming only s-wave scattering [106]:

$$\begin{eqnarray}&&{\gamma }_{t}\sim \displaystyle \frac{1}{\sqrt{6}}\sigma ({u}_{f}){n}_{f}{u}_{f},\end{eqnarray} \tag{ 99 }$$

where n_f is the number density of foreign atoms, u_f is their average velocity, and $\sigma ({u}_{f})$ is the cross-section for the process. Using kinetic theory, we can approximate these quantities by [105]:

$$\begin{eqnarray}&&{u}_{f}=\sqrt{2{k}_{B}T/{m}_{f}},\ \ {n}_{f}=P/({k}_{B}T),\ \ {\sigma }_{f}=7.57{\left(1.033\right)}^{2}{\left(\displaystyle \frac{{C}_{6}}{{\hslash }{u}_{f}}\right)}^{2/5},\end{eqnarray} \tag{ 100 }$$

where P is the pressure of the vacuum chamber and C₆ is the Van der Waals constants from the Van der Waals potential $V(r)=-{C}_{6}/{r}^{6}$. Various interaction cross-sections have been calculated for these processes in [105] assuming T is room temperature.

The N-particle correlation for a NOON state under the master equation (98) is then:

$$\begin{eqnarray}&&\langle {\hat{a}}_{L}^{\dagger N}{\hat{a}}_{R}^{N}\rangle (t)={{\rm{e}}}^{-{\gamma }_{f}{Nt}}\langle {\hat{a}}_{L}^{\dagger N}{\hat{a}}_{R}^{N}{\rangle }_{0},\end{eqnarray} \tag{ 101 }$$

so that we require

$$\begin{eqnarray}&&\displaystyle \frac{{{Gm}}^{2}N}{R{\hslash }}\gg {\gamma }_{f}.\end{eqnarray} \tag{ 102 }$$

Note that the environmental decoherence rate here is equal to the atomic loss rate. This is because, for an exact NOON state, the loss of one atom means that the density operator is now a mixture of the states $| ({N}-1)0\rangle$ and $| 0({N}-1)\rangle$, neither of which is itself a NOON state [39]. Therefore, one scattering event is enough to collapse the NOON state into all-left or all-right states. However, in practice it is unlikely that an exact NOON state will be formed, and instead a more general macroscopic superposition state, such as $| \psi \rangle =(| {N}/10,9{N}/10\rangle +| 9{N}/10,{N}/10\rangle )/\sqrt{2}$, would be more probable. Indeed, these type of states would be formed in the process where the scattering length is suddenly changed to a negative value [39]. Single-atom losses for these states would still result in similarly 'good' macroscopic superposition states, such that the effect of scattering a foreign atom may not have a significant detrimental effect [39]. Of course we also need to determine how the GQSR rate might change for the approximate NOON states, and this will be the concern of future work.

3.3.4. Decoherence from the trapping potential

Optical, magnetic or opto-magnetic traps can be used for the implementation of the double-well potential. These electromagnetic traps can also cause decoherence of a NOON state. For example, in an optical trap, decoherence of NOON states can come from spontaneously scattered/diffracted photons [39, 107], and phase noise of the laser beam [108]. However, to generate the required large numbers of atoms, it is likely that a pure magnetic trap should be used, such as that in [109, 110]. Surprisingly, decoherence of macroscopic superposition states due to fluctuations of a magnetic field has been found to be independent of the total particle number [111], improving the feasibility of generating such states.

We have investigated testing a unified theory of GR and QT with a BEC. In particular, we have considered testing a proposal for a unified theory that is based on the 'gravitizing QT' approach rather than the conventional 'quantizing gravity' approach. In section 2, we examined how, if we attempt to make QT consistent with the equivalence principle of GR, then a possible resolution is to consider making modifications to QT that would lead to a violation of the superposition principle of QT where the degree of violation is dependent on the gravitational interaction and configuration of the system. Since this increases for more massive systems, the proposal can provide an objective state reduction that is consistent with current experiments, thus resolving the measurement problem of QT, which would, on other hand, be expected to persist for the 'quantizing gravity' approach and conventional quantum gravity theories. QT is predicted to breakdown when the mass of a quantum system is near the Planck mass scale, allowing for experimental tests that are far more achievable than those generally required for distinguishing conventional quantum gravity theories, where the relevant effects are anticipated near the Planck length scale.

In the proposal considered here for a unified theory of GR and QT, quantum superposition states are expected to decay to localized states with an average lifetime that is (in the Newtonian limit) reciprocally related to the self-energy of the difference between the mass distributions of the localized states, E_G [3], which is dependent on the mass and configuration of the system. This has been generically considered for displaced, uniform, spherical mass distributions. However, BECs tend to have non-uniform mass distributions, and so we have extended this to the quadratic and Gaussian density distributions that are usually found in BEC experiments, but which may also be applicable to other systems. Since they are often generated in BEC experiments, we have also considered non-uniform spheroidal mass distributions, as well as uniform ones that would be approximated in nano/micro-object experiments, finding that the average lifetime of state reduction can be reduced compared to the spherical case. Due to the particular dependence that the GQSR considered here has on the geometry of the superposed objects, this analysis could also be used to distinguish the GQSR from other, and potentially non-gravitational, collapse models, such as the continuous spontaneous localization (CSL) model [112].

To probe the GQSR, we have considered a BEC in a double-well trap that is placed in a macroscopic superposition state of two locations. Assuming that the state reduction is a Poisson process similar to particle decay, we have compared the rate of wavefunction collapse against prominent channels of environmental decoherence in BEC systems. For the rate of decoherence to be significantly less than the rate of collapse, we estimate that the BEC should have greater than 10⁸ or 10⁹ atoms, depending on the choice of a free parameter in the GQSR proposal, and that the scattering length is reduced using an external magnetic field while maintaining a macroscopic superposition state. Being able to control the atom-atom interactions provides a unique asset to BEC tests.

We have concentrated on exact NOON states for estimating collapse and decoherence rates. However, as with experimental proposals based on nano/micro-objects, these states would be challenging to create and approximations to these states are more likely to be generated in experiments. Although estimating environmental decoherence for approximate NOON states is a relatively simple task, the GQSR needs to be extended to be able to handle these states. One possibility is to follow the approach of Diósi [15], but there may be other, more general, alternatives, which will be the concern of future work. We have also concentrated on only three-dimensional BECs, but prolate and oblate BECs with high ellipticity could move into a quasi-one and -two dimensional regime, potentially reducing environmental decoherence processes such as that from three-body recombination [97–99]. In this case, environmental decoherence could be reduced relative to spherical BECs, whereas the collapse rate would be increased, improving the feasibility of experimental tests.

If signals of this proposal were not observed in experiments, then, depending on the achievable experimental parameters, this could place severe constraints on the model (for example, the value of γ), and potentially rule it out. It would also likely place constraints on other models of objective state reduction, such as CSL [112]. However, if signals were observed, then we would have the first evidence of how GR and QT must combine to form a consistent, unified theory. Furthermore, it would explain the mysterious measurement process in QT and provide it with a well-defined classical limit.

We thank participants of the Gravity in the lab 2018 workshop in Benasque, Spain for useful discussions and comments. In particular, we would like to thank Philippe Bouyer, Devang Naik, Hendrik Ulbricht, Marko Toroš Daniel Goldwater and Michael Tolan. RP would like to thank John Moussouris for a personal endowment through the University of Oxford. RH and IF would like to acknowledge the support of the Penrose Institute.

Here we calculate the self-energy of the difference between uniform, displaced spherical mass distributions. Taking the radius of the sphere states to be R and their mass M, then their density functions are defined by

$$\begin{eqnarray}\rho ({\boldsymbol{r}})=\left\{\begin{array}{cl}{\rho }_{0}:= M/\left(\displaystyle \frac{4}{3}\pi {R}^{3}\right) & \mathrm{if}\ {\boldsymbol{r}}\mathrm{is}\ \mathrm{inside}\ \mathrm{sphere},\\ 0 & \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ A.1 }$$

In terms of the step function θ(x), we can write the density functions as

$$\begin{eqnarray}&&\rho ({\boldsymbol{r}})={\rho }_{0}\theta (R-{r}_{\rho })\theta \left({R}^{2}-{r}_{\rho }^{2}-{z}^{2}\right),\end{eqnarray} \tag{ A.2 }$$

$$\begin{eqnarray}&&\rho ^{\prime} ({\boldsymbol{r}})={\rho }_{0}\theta (R-{r}_{\rho })\theta \left({R}^{2}-{r}_{\rho }^{2}-{\left(z-b\right)}^{2}\right),\end{eqnarray} \tag{ A.3 }$$

where r_ρ is the cylindrical radial coordinate, and we have taken the sphere states to be displaced along the z coordinate by a distance b, with the $\rho ({\boldsymbol{r}})$ sphere state being at the origin of our coordinate system.

These density functions can be plugged into (19). However, we find it simpler to work with (18) to calculate E_G for this situation. In this case we need the gravitational potential inside and outside a sphere:

$$\begin{eqnarray}\phi =\left\{\begin{array}{ll}{\phi }_{{\rm{in}}}:= -\displaystyle \frac{{GM}}{R}\left(\displaystyle \frac{3}{2}-\displaystyle \frac{{r}^{2}}{2{R}^{2}}\right) & \mathrm{if}\ r\leqslant R,\\ {\phi }_{{\rm{out}}}:= -\displaystyle \frac{{GM}}{r} & \mathrm{if}\ r\geqslant R,\end{array}\right.\end{eqnarray} \tag{ A.4 }$$

where now we are working in spherical coordinates (r, θ, ψ). Taking γ = 1/(8π) then, due to the symmetry of the problem, we can use

$$\begin{eqnarray}&&{E}_{G}=\displaystyle \frac{1}{2}\int (\phi -\phi ^{\prime} )(\rho ^{\prime} -\rho ){{\rm{d}}}^{3}{\boldsymbol{r}}\end{eqnarray} \tag{ A.5 }$$

$$\begin{eqnarray}&&=\int \phi (\rho ^{\prime} -\rho ){{\rm{d}}}^{3}{\boldsymbol{r}}.\end{eqnarray} \tag{ A.6 }$$

We first consider the term $\phi \rho ^{\prime}$, which is related to the gravitational interaction energy [3]. When b > 2R then, following Gauss's law, the gravitational interaction energy is simply −GM²/b. We can calculate this by choosing the origin of our coordinate system to be at the center of the $\rho ({\boldsymbol{r}})$ sphere and integrate its potential over the density of the other sphere using surfaces of constant radial coordinate r:

$$\begin{eqnarray}&&\int \phi \rho ^{\prime} {{\rm{d}}}^{3}{\boldsymbol{r}}=-{\rho }_{0}{\int }_{0}^{2\pi }{\int }_{0}^{{\cos }^{-1}({z}_{r}/r)}{\int }_{b-R}^{b+R}{\phi }_{{out}}{r}^{2}\sin \theta {\rm{d}}{r}{\rm{d}}\theta {\rm{d}}\psi =-\displaystyle \frac{{{GM}}^{2}}{b}.\end{eqnarray} \tag{ A.7 }$$

Using the same method for $R\leqslant b\leqslant 2R$, we find:

$$\begin{eqnarray}&&\int \phi \rho ^{\prime} {{\rm{d}}}^{3}{\boldsymbol{r}}=-{\rho }_{0}{\int }_{0}^{2\pi }{\int }_{0}^{{\cos }^{-1}({z}_{r}/r)}\left[{\int }_{b-R}^{R}{\phi }_{{\rm{in}}}+{\int }_{R}^{b+R}{\phi }_{{\rm{out}}}\right]{r}^{2}\sin \theta {\rm{d}}{r}{\rm{d}}\theta {\rm{d}}\psi \end{eqnarray} \tag{ A.8 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{{{GM}}^{2}}{R}\left(-\displaystyle \frac{6}{5}+2{\lambda }^{2}-\displaystyle \frac{3}{2}{\lambda }^{3}+\displaystyle \frac{1}{5}{\lambda }^{5}\right).\end{eqnarray} \tag{ A.9 }$$

Finally, for $b\leqslant R$, we have:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \int \phi \rho ^{\prime} {{\rm{d}}}^{3}{\boldsymbol{r}} & = & -{\rho }_{0}{\displaystyle \int }_{0}^{2\pi }\left[\left({\displaystyle \int }_{0}^{\pi }{\displaystyle \int }_{0}^{R-b}+{\displaystyle \int }_{0}^{{\cos }^{-1}({z}_{r}/r)}{\displaystyle \int }_{R-b}^{R}\right){\phi }_{{\rm{in}}}\right.\\ & & \left.+{\displaystyle \int }_{0}^{{\cos }^{-1}({z}_{r}/r)}{\displaystyle \int }_{R+b}^{R}{\phi }_{{\rm{out}}}\right]{r}^{2}\sin \theta {\rm{d}}{r}{\rm{d}}\theta {\rm{d}}\psi \end{array}\end{eqnarray} \tag{ A.10 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{{{GM}}^{2}}{R}\left(-\displaystyle \frac{6}{5}+2{\lambda }^{2}-\displaystyle \frac{3}{2}{\lambda }^{3}+\displaystyle \frac{1}{5}{\lambda }^{5}\right),\end{eqnarray} \tag{ A.11 }$$

which, as expected is the same as the previous result. Now we consider the term ϕ ρ in (A.6). We can simply extract this from the above result when b = 0 such that:

$$\begin{eqnarray}&&\int \phi \rho {{\rm{d}}}^{3}{\boldsymbol{r}}=-\displaystyle \frac{6{{GM}}^{2}}{5R},\end{eqnarray} \tag{ A.12 }$$

which is simply twice the gravitational self-energy of a sphere. Putting this altogether, we then obtain:

$$\begin{eqnarray}{E}_{G}=\left\{\begin{array}{cl}\displaystyle \frac{{{GM}}^{2}}{R}\left(2{\lambda }^{2}-\displaystyle \frac{3}{2}{\lambda }^{3}+\displaystyle \frac{1}{5}{\lambda }^{5}\right) & \mathrm{if}\ 0\leqslant \lambda \leqslant 1,\\ \displaystyle \frac{{{GM}}^{2}}{R}\left(\displaystyle \frac{6}{5}-\displaystyle \frac{1}{2\lambda }\right) & \mathrm{if}\ \lambda \geqslant 1.\end{array}\right.\end{eqnarray} \tag{ A.13 }$$

Here we calculate the self-energy of the difference between displaced spherical BEC mass distributions in the TF approximation. Taking the radius of the sphere states to be R and their total mass M, then their density functions are defined by

$$\begin{eqnarray}&&\rho ({\boldsymbol{r}})={\rho }_{0}({\boldsymbol{r}})\theta (R-{r}_{\rho })\theta \left({R}^{2}-{r}_{\rho }^{2}-{z}^{2}\right),\end{eqnarray} \tag{ B.1 }$$

$$\begin{eqnarray}&&\rho ^{\prime} ({\boldsymbol{r}})=\rho {{\prime} }_{0}({\boldsymbol{r}})\theta (R-{r}_{\rho })\theta \left({R}^{2}-{r}_{\rho }^{2}-{\left(z-b\right)}^{2}\right),\end{eqnarray} \tag{ B.2 }$$

where

$$\begin{eqnarray}&&{\rho }_{0}({\boldsymbol{r}}):= \displaystyle \frac{5}{2}{\rho }_{0}\left(1-\left({r}_{\rho }^{2}+{z}^{2}\right)/{R}^{2}\right)=\displaystyle \frac{5}{2}{\rho }_{0}\left(1-{r}^{2}/{R}^{2}\right),\end{eqnarray} \tag{ B.3 }$$

$$\begin{eqnarray}&&\rho {{\prime} }_{0}({\boldsymbol{r}}):= \displaystyle \frac{5}{2}{\rho }_{0}\left(1-\left({r}_{\rho }^{2}+{\left[z-b\right]}^{2}\right)/{R}^{2}\right)=\displaystyle \frac{5}{2}{\rho }_{0}\left(1-r{{\prime} }^{2}/{R}^{2}\right),\end{eqnarray} \tag{ B.4 }$$

with $r{{\prime} }^{2}={r}^{2}+{b}^{2}-2{rb}\cos \theta$ and ${\rho }_{0}:= M/((4/3)\pi {R}^{3})$ as before. We again use (18) to calculate E_G for this situation. In this case we need the gravitational potential inside and outside the sphere. From Gauss's law, the outside potential is of course the same as in the uniform situation. To find the inner potential, we can also apply Gauss's law:

$$\begin{eqnarray}&&\oint {\boldsymbol{g}}.{\rm{d}}{\boldsymbol{S}}=-4\pi {{GM}}_{r},\end{eqnarray} \tag{ B.5 }$$

where we choose a spherical surface of constant radius r within the sphere such that M_r is the total mass within this spherical surface and is given by

$$\begin{eqnarray}&&{M}_{r}=4\pi {\int }_{0}^{r}\rho ({\boldsymbol{r}}^{\prime} ){r}^{{\prime} 2}{\rm{d}}{r}^{\prime} ,\end{eqnarray} \tag{ B.6 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{M}{2{R}^{5}}\left(5{R}^{2}{r}^{3}-3{r}^{5}\right),\end{eqnarray} \tag{ B.7 }$$

where M is the total mass of the spherical BEC. Therefore, the field inside the sphere is given by:

$$\begin{eqnarray}&&{\boldsymbol{g}}=-\displaystyle \frac{{GM}}{2{R}^{5}{r}^{2}}\left(5{R}^{2}{r}^{3}-3{r}^{5}\right),\end{eqnarray} \tag{ B.8 }$$

and the potential then can be found through:

$$\begin{eqnarray}&&{\phi }_{{\rm{in}}}(r)=-{\int }_{\infty }^{R}\displaystyle \frac{-{GM}}{{r}^{2}}{\rm{d}}{r}-{\int }_{R}^{r}\displaystyle \frac{-{GM}}{2{R}^{5}{r}^{2}}\left(5{R}^{2}{r}^{3}-3{r}^{5}\right)\end{eqnarray} \tag{ B.9 }$$

$$\begin{eqnarray}&&=-\displaystyle \frac{{GM}}{8{R}^{5}}\left(15{R}^{4}-10{R}^{2}{r}^{2}+3{r}^{4}\right).\end{eqnarray} \tag{ B.10 }$$

The gravitational potential of a spherical BEC in the TF approximation is then

$$\begin{eqnarray}\phi =\left\{\begin{array}{ll}{\phi }_{{\rm{in}}}:= -\displaystyle \frac{{GM}}{8R}\left(15-10{r}^{2}/{R}^{2}+3{r}^{4}/{R}^{4}\right) & \mathrm{if}\ r\leqslant R,\\ {\phi }_{{\rm{out}}}:= -\displaystyle \frac{{GM}}{r} & \mathrm{if}\ r\geqslant R.\end{array}\right.\end{eqnarray} \tag{ B.11 }$$

The rest of the calculation now proceeds similar to the uniform case. We first consider the term $\phi \rho ^{\prime}$ in (A.6). When b > 2R we find, due to Gauss's law again, that this is simply −GM²/b as before. For $0\leqslant b\leqslant 2R$ we choose the origin of our coordinate system to be at the center of the $\rho ({\boldsymbol{r}})$ sphere state and again integrate its potential over the density of the other sphere state using surfaces of constant radial coordinate r:

$$\begin{eqnarray}&&\int \phi \rho ^{\prime} {{\rm{d}}}^{3}{\boldsymbol{r}}=-{\int }_{0}^{2\pi }{\int }_{0}^{{\cos }^{-1}({z}_{r}/r)}\left[{\int }_{b-R}^{R}{\phi }_{{\rm{in}}}+{\int }_{R}^{b+R}{\phi }_{{\rm{out}}}\right]\rho {{\prime} }_{0}({\boldsymbol{r}}){r}^{2}\sin \theta {\rm{d}}{r}{\rm{d}}\theta {\rm{d}}\psi \end{eqnarray} \tag{ B.12 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{{{GM}}^{2}}{R}\left(-\displaystyle \frac{10}{7}+\displaystyle \frac{20}{7}{\lambda }^{2}-6{\lambda }^{4}+5{\lambda }^{6}-\displaystyle \frac{15}{14}{\lambda }^{7}+\displaystyle \frac{1}{7}{\lambda }^{9}\right).\end{eqnarray} \tag{ B.13 }$$

For the term ϕ ρ in (A.6) we can again simply extract this from the above result when b = 0:

$$\begin{eqnarray}&&\int \phi \rho =-\displaystyle \frac{10{{GM}}^{2}}{7R}.\end{eqnarray} \tag{ B.14 }$$

Putting this altogether, we obtain:

$$\begin{eqnarray}{E}_{G}=\left\{\begin{array}{ll}\displaystyle \frac{{{GM}}^{2}}{R}\left(\displaystyle \frac{20}{7}{\lambda }^{2}-6{\lambda }^{4}+5{\lambda }^{5}-\displaystyle \frac{15}{14}{\lambda }^{7}+\displaystyle \frac{1}{7}{\lambda }^{9}\right) & \mathrm{if}\ 0\leqslant \lambda \leqslant 1,\\ \displaystyle \frac{{{GM}}^{2}}{R}\left(\displaystyle \frac{10}{7}-\displaystyle \frac{1}{2\lambda }\right) & \mathrm{if}\ \lambda \geqslant 1.\end{array}\right.\end{eqnarray} \tag{ B.15 }$$

Here we calculate the self-energy of the difference between displaced spherical BEC mass distributions in the Gaussian approximation. Taking the sphere to have total mass M, then the density functions are defined by

$$\begin{eqnarray}&&\rho ({\boldsymbol{r}})=\displaystyle \frac{4}{3\sqrt{\pi }}{\rho }_{0}{{\rm{e}}}^{-{r}^{2}/{R}_{0}^{{\prime} 2}},\end{eqnarray} \tag{ C.1 }$$

$$\begin{eqnarray}&&\rho ^{\prime} ({\boldsymbol{r}})=\displaystyle \frac{4}{3\sqrt{\pi }}{\rho }_{0}{{\rm{e}}}^{-r{{\prime} }^{2}/{R}_{0}^{{\prime} 2}},\end{eqnarray} \tag{ C.2 }$$

where ${\rho }_{0}:= M/((4/3)\pi {R}_{0}^{{\prime} 3})$ and ${R}_{0}^{{\prime} }$ is given by (41) and can be taken as a measure for the size of the condensate [66]. However, we do not take this as a discontinuous cut in the density and instead keep the the wavefunction of the condensate has having infinite extent. Following the procedure outlined in the previous section using Gauss's law (or using the method outlined in appendix F), the potential of a Gaussian sphere is found to be

$$\begin{eqnarray}&&\phi (r)=-\displaystyle \frac{{GM}}{r}\mathrm{erf}(r/{R}_{0}^{{\prime} }),\end{eqnarray} \tag{ C.3 }$$

where erf(x) is the error function. This tends to $-2{GM}/\left(\sqrt{\pi }{R}_{0}^{{\prime} }\right)$ in the limit that $r\to 0$. The ϕ ρ term of (A.6) is then found to be

$$\begin{eqnarray}&&\int {{\rm{d}}}^{3}{\boldsymbol{r}}\phi \rho ={\int }_{0}^{2\pi }{\int }_{0}^{\pi }{\int }_{0}^{\infty }\rho ({\boldsymbol{r}})\phi (r)=-\displaystyle \frac{{{GM}}^{2}}{{R}_{0}^{{\prime} }}\sqrt{\displaystyle \frac{2}{\pi }},\end{eqnarray} \tag{ C.4 }$$

and the $\phi \rho ^{\prime}$ term is

$$\begin{eqnarray}&&\int {{\rm{d}}}^{3}{\boldsymbol{r}}\phi \rho ^{\prime} ={\int }_{0}^{2\pi }{\int }_{0}^{\pi }{\int }_{0}^{\infty }\rho ^{\prime} ({\boldsymbol{r}})\phi (r)=-\displaystyle \frac{{{GM}}^{2}}{b}\mathrm{erf}\left(\displaystyle \frac{b}{\sqrt{2}{R}_{0}^{{\prime} }}\right),\end{eqnarray} \tag{ C.5 }$$

where we have used the identity [113]:

$$\begin{eqnarray}&&{\int }_{0}^{\infty }\left[{{\rm{e}}}^{-{\left(x-A\right)}^{2}}-{{\rm{e}}}^{-{\left(x+A\right)}^{2}}\right]\mathrm{erf}(x){\rm{d}}{x}\equiv \sqrt{\pi }\mathrm{erf}\left(A/\sqrt{2}\right).\end{eqnarray} \tag{ C.6 }$$

Putting this altogether, we obtain:

$$\begin{eqnarray}&&{E}_{G}=\displaystyle \frac{{{GM}}^{2}}{{R}_{0}^{{\prime} }}\left(\sqrt{\displaystyle \frac{2}{\pi }}-\displaystyle \frac{1}{2\lambda }\mathrm{erf}\left(\sqrt{2}\lambda \right)\right)\end{eqnarray} \tag{ C.7 }$$

where we have defined $\lambda := b/(2{R}_{0}^{{\prime} })$.

Here we consider the self-energy of the difference between uniform, displaced spheroidal mass distributions. Following the previous sections, we work with (18) to calculate E_G. In this case we need the gravitational potential inside and outside the spheroid. This is simplest in spheroidal coordinates: in prolate spheroidal coordinates, the gravitational potential of prolate spheroid is [114, 115]:

$$\begin{eqnarray}\phi =\left\{\begin{array}{ll}{\phi }_{{\rm{in}}}^{{\rm{prolate}}}:= -\displaystyle \frac{{GM}}{l}\left({Q}_{0}({\xi }_{0}){G}_{1}(\xi ,\eta )+\displaystyle \frac{\xi \left({\xi }^{2}-1\right)}{{\xi }_{0}\left[{\xi }_{0}^{2}-1\right]}{G}_{2}\left(\xi ,\eta \right)+\displaystyle \frac{{\xi }_{0}^{2}-{\xi }^{2}}{{\xi }_{0}\left[{\xi }_{0}^{2}-1\right]}{G}_{3}\left(\xi ,\eta \right)\right) & \mathrm{if}\ \xi \leqslant {\xi }_{0}\\ {\phi }_{{\rm{out}}}^{{\rm{prolate}}}:= -\displaystyle \frac{{GM}}{l}\left({Q}_{0}(\xi )-{Q}_{2}(\xi ){P}_{2}(\eta )\right) & \mathrm{if}\xi \geqslant {\xi }_{0},\end{array}\right.\end{eqnarray} \tag{ D.1 }$$

where

$$\begin{eqnarray}&&{G}_{1}(\xi ,\eta ):= 1-{P}_{2}(\xi ){P}_{2}(\eta ),\end{eqnarray} \tag{ D.2 }$$

$$\begin{eqnarray}&&{G}_{2}(\xi ,\eta ):= \displaystyle \frac{3}{2}{P}_{2}(\eta )\xi ,\end{eqnarray} \tag{ D.3 }$$

$$\begin{eqnarray}&&{G}_{3}(\xi ,\eta ):= \displaystyle \frac{1}{2}+{P}_{2}(\xi ){P}_{2}(\eta ),\end{eqnarray} \tag{ D.4 }$$

P_n(x) and Q_n(x) are Legendre polynomials of the first and second kind¹⁷ :

$$\begin{eqnarray}&&{P}_{2}(x)=\displaystyle \frac{1}{2}(3{x}^{2}-1),\end{eqnarray} \tag{ D.6 }$$

$$\begin{eqnarray}&&{Q}_{0}(x)=(1/2)\mathrm{ln}\left(\displaystyle \frac{x+1}{x-1}\right)={\tan }^{-1}x+\displaystyle \frac{1}{2}{\rm{i}}\pi ,\end{eqnarray} \tag{ D.7 }$$

$$\begin{eqnarray}&&{Q}_{2}(x)={P}_{2}(x){Q}_{0}(x)-\displaystyle \frac{3}{2}x;\end{eqnarray} \tag{ D.8 }$$

$l=\sqrt{{c}^{2}-{a}^{2}}$ is the focal distance with a and c the equatorial and polar radii respectively (which are respectively the semi-minor and semi-major axes for the prolate but semi-major and semi-minor for the oblate case); (ξ, η, ψ) are prolate spheroidal coordinates with

$$\begin{eqnarray}&&x/l=\sqrt{{\xi }^{2}-1}\sin \nu \cos \psi ,\end{eqnarray} \tag{ D.9 }$$

$$\begin{eqnarray}&&y/l=\sqrt{{\xi }^{2}-1}\sin \nu \sin \psi ,\end{eqnarray} \tag{ D.10 }$$

$$\begin{eqnarray}&&z/l=\xi \cos \nu ,\end{eqnarray} \tag{ D.11 }$$

using $\eta =\cos \nu ;$ and ξ₀ := c/l is the value of ξ at the surface of the prolate spheroid. For the potential of the oblate spheroid, just replace ξ with iξ, l by −il and ξ₀ with iξ₀ [114, 115].

We first consider the term ϕρ in (18) for a prolate spheroid:

$$\begin{eqnarray}&&\int {{\rm{d}}}^{3}{\boldsymbol{r}}\phi \rho ={\rho }_{0}^{{\rm{spheroid}}}{\int }_{0}^{2\pi }{\int }_{-1}^{1}{\int }_{1}^{{\xi }_{0}}{l}^{3}({\xi }^{2}-{\eta }^{2}){\phi }_{{\rm{in}}}^{{\rm{prolate}}}{\rm{d}}\xi {\rm{d}}\eta {\rm{d}}\psi \end{eqnarray} \tag{ D.12 }$$

$$\begin{eqnarray}&&=\displaystyle \frac{6{{GM}}^{2}}{5l}{\coth }^{-1}{\xi }_{0},\end{eqnarray} \tag{ D.13 }$$

$$\begin{eqnarray}&&=\displaystyle \frac{6{{GM}}^{2}}{5l}{\tanh }^{-1}e,\end{eqnarray} \tag{ D.14 }$$

where e is the ellipticity (e := c/l for the prolate case) and ${\rho }_{0}^{{\rm{spheroid}}}$ = M/((4/3)π${a}^{2}c)$ is the density of a uniform spheroid. For an oblate spheroid, we just need to replace ξ₀ with iξ₀ to obtain

$$\begin{eqnarray}&&\displaystyle \frac{6{{GM}}^{2}}{5l}{\cot }^{-1}{\xi }_{0}=\displaystyle \frac{6{{GM}}^{2}}{5l}{\sin }^{-1}e,\end{eqnarray} \tag{ D.15 }$$

where e := a/l is now the ellipticity of an oblate spheroid with $l=\sqrt{{a}^{2}-{c}^{2}}$ its focal distance.

We now consider the term $\phi \rho ^{\prime}$ in (18). For the prolate spheroid in configuration (b) in figure 3, we choose to integrate over surfaces of constant ξ. When $b\geqslant 2c$, we use

$$\begin{eqnarray}&&\int \phi \rho ^{\prime} {{\rm{d}}}^{3}{\boldsymbol{r}}={\rho }_{0}^{{\rm{spheroid}}}{\int }_{0}^{2\pi }{\int }_{{\eta }_{\xi }}^{1}{\int }_{(b-c)/l}^{(b+c)/l}{l}^{3}({\xi }^{2}-{\eta }^{2}){\phi }_{{\rm{out}}}^{{\rm{prolate}}}{\rm{d}}\xi {\rm{d}}\eta {\rm{d}}\psi ,\end{eqnarray} \tag{ D.16 }$$

where η_ξ is the η-coordinate where the ξ-surface meets with the surface of the $\rho ^{\prime}$ spheroid. When $0\leqslant b\leqslant 2c$, we use

$$\begin{eqnarray}&&\displaystyle \int \phi \rho ^{\prime} {{\rm{d}}}^{3}{\boldsymbol{r}}={\rho }_{0}^{{\rm{spheroid}}}{\displaystyle \int }_{0}^{2\pi }{\displaystyle \int }_{{\eta }_{\xi }}^{1}\left[{\displaystyle \int }_{(b-c)/l}^{c/l}{\phi }_{{\rm{in}}}^{{\rm{prolate}}}\right.\end{eqnarray} \tag{ D.17 }$$

$$\begin{eqnarray}&&\left.+\ {\displaystyle \int }_{c/l}^{(b+c)/l}{\phi }_{{\rm{out}}}^{{\rm{prolate}}}\right]{l}^{3}({\xi }^{2}-{\eta }^{2}){\rm{d}}\xi {\rm{d}}\eta {\rm{d}}\psi .\end{eqnarray} \tag{ D.18 }$$

The result for E_G to first order in , where a = c, for a prolate spheroid is provided in section 2.2.

Another option, which is more suited to an oblate spheroidal coordinate system, is to integrate over surfaces of constant η. Choosing now to work instead with an oblate spheroid coordinate system then, when $b\geqslant 2c$, we use

$$\begin{eqnarray}&&\int \phi \rho ^{\prime} {{\rm{d}}}^{3}{\boldsymbol{r}}={\rho }_{0}^{{\rm{spheroid}}}{\int }_{0}^{2\pi }{\int }_{{\eta }_{{\rm{\max }}}}^{1}{\int }_{{\xi }_{1}}^{{\xi }_{2}}{l}^{3}({\xi }^{2}+{\eta }^{2}){\phi }_{{\rm{out}}}^{{\rm{oblate}}}{\rm{d}}\xi {\rm{d}}\eta {\rm{d}}\psi ,\end{eqnarray} \tag{ D.19 }$$

where ξ₁ and ξ₂ are the two values of ξ where the constant η-surface crosses the $\rho ^{\prime}$ spheroid state, and η_max is the value of η where there is only one ξ solution i.e. ξ₁ = ξ₂. When $0\leqslant b\leqslant 2c$, we use

$$\begin{eqnarray}&&\int \phi \rho ^{\prime} {{\rm{d}}}^{3}{\boldsymbol{r}}={\rho }_{0}^{{\rm{spheroid}}}{\int }_{0}^{2\pi }\left[{\int }_{{\eta }_{{\rm{int}}}}^{1}\left({\int }_{{\xi }_{1}}^{{\xi }_{0}}{\phi }_{{\rm{in}}}^{{\rm{oblate}}}+{\int }_{{\xi }_{0}}^{{\xi }_{2}}{\phi }_{{\rm{out}}}^{{\rm{oblate}}}\right)\right.\end{eqnarray} \tag{ D.20 }$$

$$\begin{eqnarray}&&\left.+\ {\int }_{{\eta }_{{\rm{\max }}}}^{{\eta }_{{\rm{int}}}}{\int }_{{\xi }_{1}}^{{\xi }_{2}}{\phi }_{{\rm{out}}}^{{\rm{oblate}}}\right]{l}^{3}\left({\xi }^{2}+{\eta }^{2}\right){\rm{d}}\xi {\rm{d}}\eta {\rm{d}}\psi ,\end{eqnarray} \tag{ D.21 }$$

where η_int = b/(2l ξ₀) is the η-coordinate where the two spheroids meet. Once the result for E_G is obtained, the prolate spheroid case can be found via ${\xi }_{0}\to {\rm{i}}{\xi }_{0}$ as above. When $\epsilon \ll 1$, where now c := a, E_G for an oblate spheroid displaced along its symmetry axis can be approximated by (25) in section 2.2.

Unlike in the spherical case, equipotential surfaces are not similar-shaped spheroids (or confocal spheroids), emphasizing that Gauss's law is not as useful for these objects. Therefore, integrating over constant ξ or η surfaces is not as simple. An alternative is to use cylindrical coordinates where the prolate spheroid potential inside and outside the spheroid is given by:

$$\begin{eqnarray}&&{\phi }_{{\rm{in}}}^{{\rm{prolate}}}:= -\displaystyle \frac{3{GM}}{4{l}^{3}}\left[\left(2{l}^{2}+D\right){\mathrm{csch}}^{-1}\left(\displaystyle \frac{a}{l}\right)-\displaystyle \frac{l\left({c}^{2}{r}^{2}-2{a}^{2}{z}^{2}\right)}{{a}^{2}c}\right],\end{eqnarray} \tag{ D.22 }$$

$$\begin{eqnarray}&&{\phi }_{{\rm{out}}}^{{\rm{prolate}}}:= -\displaystyle \frac{3{GM}}{4{l}^{3}}\left[\left(2{l}^{2}+D\right){\sinh }^{-1}\left(\displaystyle \frac{\sqrt{2}l}{{B}_{p}}\right)-\displaystyle \frac{\sqrt{2}l\left({A}_{p}D+{l}^{2}C\right)}{{E}_{p}{B}_{p}^{2}}\right],\end{eqnarray} \tag{ D.23 }$$

and the oblate potential is found by taking $l\to -{\rm{i}}{l}$, to obtain:

$$\begin{eqnarray}&&{\phi }_{{\rm{in}}}^{{\rm{oblate}}}:= -\displaystyle \frac{3{GM}}{4{l}^{3}}\left[\left(2{l}^{2}-D\right){\csc }^{-1}\left(\displaystyle \frac{a}{l}\right)+\displaystyle \frac{l\left({c}^{2}{r}^{2}-2{a}^{2}{z}^{2}\right)}{{a}^{2}c}\right],\end{eqnarray} \tag{ D.24 }$$

$$\begin{eqnarray}&&{\phi }_{{\rm{out}}}^{{\rm{oblate}}}:= -\displaystyle \frac{3{GM}}{4{l}^{3}}\left[\left(2{l}^{2}-D\right){\sin }^{-1}\left(\displaystyle \frac{\sqrt{2}l}{{B}_{o}}\right)+\displaystyle \frac{\sqrt{2}l\left({A}_{o}D-{l}^{2}C\right)}{{E}_{o}{B}_{o}^{2}}\right],\end{eqnarray} \tag{ D.25 }$$

where

$$\begin{eqnarray}&&{A}_{p}:= {r}^{2}+{z}^{2}+\sqrt{{l}^{4}+2{l}^{2}\left({r}^{2}-{z}^{2}\right)+{\left({z}^{2}+{r}^{2}\right)}^{2}},\end{eqnarray} \tag{ D.26 }$$

$$\begin{eqnarray}&&{A}_{o}:= {r}^{2}+{z}^{2}+\sqrt{{z}^{4}+2{z}^{2}\left({r}^{2}+{l}^{2}\right)+{\left({l}^{2}-{r}^{2}\right)}^{2}},\end{eqnarray} \tag{ D.27 }$$

$$\begin{eqnarray}&&{B}_{p}:= \sqrt{{A}_{p}-{l}^{2}},\,{E}_{p}:= \sqrt{{A}_{p}+{l}^{2}},\end{eqnarray} \tag{ D.28 }$$

$$\begin{eqnarray}&&{B}_{o}:= \sqrt{{A}_{o}+{l}^{2}},\,{E}_{o}:= \sqrt{{A}_{o}-{l}^{2}},\end{eqnarray} \tag{ D.29 }$$

$$\begin{eqnarray}&&C:= {r}^{2}+2{z}^{2},\end{eqnarray} \tag{ D.30 }$$

$$\begin{eqnarray}&&D:= {r}^{2}-2{z}^{2}.\end{eqnarray} \tag{ D.31 }$$

To second order in ellipticity, the potentials, in spherical coordinates, become:

$$\begin{eqnarray}\phi =\left\{\begin{array}{cl}{\phi }_{{\rm{in}}}:= {\phi }_{{\rm{in}}}^{{\rm{sphere}}}\pm \displaystyle \frac{{GM}}{5{R}^{3}}{P}_{2}(\cos \theta ){e}^{2} & \mathrm{if}\ r\leqslant r(\theta ),\\ {\phi }_{{\rm{out}}}:= {\phi }_{{\rm{out}}}^{{\rm{sphere}}}\pm \displaystyle \frac{{{GMR}}^{2}}{5{r}^{3}}{P}_{2}(\cos \theta ){e}^{2} & \mathrm{if}\ r\geqslant r(\theta ),\end{array}\right.\end{eqnarray} \tag{ D.32 }$$

where + is for the prolate case, − is for the oblate case, ${\phi }_{{\rm{in}}/{\rm{out}}}^{{\rm{sphere}}}$ is given by (A.4), $r(\theta ):= c{\left[1-{e}^{2}{\cos }^{2}\theta \right]}^{-1/2}$ (with the respective definitions of ellipticity for the two spheroidal cases), and we have taken both spheroids to have the same volume as a sphere with radius R.

Using the full potentials in cylindrical coordinates, for the spheroidal cases in configuration (a) and (b) in figure 3, the ϕ ρ term is then:

$$\begin{eqnarray}&&\int \phi \rho {{\rm{d}}}^{3}{\boldsymbol{r}}={\rho }_{0}^{{\rm{spheroid}}}{\int }_{0}^{2\pi }{\int }_{-c}^{c}{\int }_{0}^{a\sqrt{1-{z}^{2}/{c}^{2}}}{\phi }_{{\rm{in}}}^{{\rm{spheroid}}}{r}{\rm{d}}{r}{\rm{d}}{z}{\rm{d}}\psi .\end{eqnarray} \tag{ D.33 }$$

where ${\phi }_{{\rm{in}}}^{{\rm{spheroid}}}$ is given by either (D.22) or (D.24). For $b\geqslant 2c$, the $\phi \rho ^{\prime}$ term is:

$$\begin{eqnarray}&&\int \phi \rho ^{\prime} {{\rm{d}}}^{3}{\boldsymbol{r}}={\rho }_{0}^{{\rm{spheroid}}}{\int }_{0}^{2\pi }{\int }_{b-c}^{b+c}{\int }_{0}^{a\sqrt{1-{\left(z-b\right)}^{2}/{c}^{2}}}{\phi }_{{\rm{out}}}^{{\rm{spheroid}}}{r}{\rm{d}}{r}{\rm{d}}{z}{\rm{d}}\psi ,\end{eqnarray} \tag{ D.34 }$$

where ${\phi }_{{\rm{out}}}^{{\rm{spheroid}}}$ is given by either (D.23) or (D.25); and for $0\leqslant b\leqslant 2c$, we can use

$$\begin{eqnarray}&&{\rho }_{0}^{{\rm{spheroid}}}{\int }_{0}^{2\pi }\left[\left({\int }_{b-c}^{b/2}{\int }_{0}^{a\sqrt{1-{\left(z-b\right)}^{2}/{c}^{2}}}+{\int }_{b/2}^{c}{\int }_{0}^{a\sqrt{1-{z}^{2}/{c}^{2}}}\right){\phi }_{{\rm{in}}}^{{\rm{spheroid}}}\right.\end{eqnarray} \tag{ D.35 }$$

$$\begin{eqnarray}&&\left.+\ \left({\int }_{b/2}^{c}{\int }_{a\sqrt{1-{z}^{2}/{c}^{2}}}^{a\sqrt{1-{\left(z-b\right)}^{2}/{c}^{2}}}+{\int }_{c}^{b+c}{\int }_{0}^{a\sqrt{1-{\left(z-b\right)}^{2}/{c}^{2}}}\right){\phi }_{{\rm{out}}}^{{\rm{spheroid}}}\right]{r}{\rm{d}}{r}{\rm{d}}{z}{\rm{d}}\psi .\end{eqnarray} \tag{ D.36 }$$

For the oblate and prolate configurations (c) and (d) in figure 3, the above procedure is just slightly modified.

Within the TF approximation, the density of spheroidal BECs is given by (see (48)):

$$\begin{eqnarray}&&{\rho }_{0}^{{\rm{oblate}}}({\boldsymbol{r}})=\displaystyle \frac{5}{2}{\rho }_{0}^{{\rm{spheroid}}}(1-{r}_{\rho }^{2}/{a}^{2}-{z}^{2}/{c}^{2})\end{eqnarray} \tag{ E.1 }$$

$$\begin{eqnarray}&&=\displaystyle \frac{5}{2}{\rho }_{0}^{{\rm{spheroid}}}\displaystyle \frac{{\xi }_{0}^{2}-{\xi }^{2}}{{\xi }_{0}^{2}[{\xi }_{0}^{2}+1]}[{\xi }_{0}^{2}+{\eta }^{2}]\end{eqnarray} \tag{ E.2 }$$

$$\begin{eqnarray}&&=\displaystyle \frac{5{\epsilon }^{2}}{2}{\rho }_{0}^{{\rm{spheroid}}}\left(1-e{{\prime} }^{2}\left[{\xi }^{2}+{\eta }^{2}\left(1-e{{\prime} }^{2}{\xi }^{2}\right)\right]\right),\end{eqnarray} \tag{ E.3 }$$

$$\begin{eqnarray}&&{\rho }_{0}^{{\rm{prolate}}}({\boldsymbol{r}})=\displaystyle \frac{5}{2{\epsilon }^{2}}{\rho }_{0}^{{\rm{spheroid}}}\left(1-{e}^{2}\left[{\xi }^{2}+{\eta }^{2}\left(1-{e}^{2}{\xi }^{2}\right)\right]\right)\end{eqnarray} \tag{ E.4 }$$

for the respective coordinate systems, where $e^{\prime} := l/c$ is the second ellipticity for an oblate spheroid.

We now find the gravitational potential of these spheroidal BECs by summing the individual potentials from point-like sources of mass ${\rm{d}}{m}=\rho ({\boldsymbol{r}}^{\prime} ){{\rm{d}}}^{3}{\boldsymbol{r}}^{\prime}$ where ${{\rm{d}}}^{3}{\boldsymbol{r}}^{\prime}$ is the volume element of the spheroid, and $\rho ({\boldsymbol{r}})$ is its density function, i.e. we use (16):

$$\begin{eqnarray}&&\phi ({\boldsymbol{x}})\equiv -G\int \displaystyle \frac{\rho ({\boldsymbol{r}}^{\prime} )}{r}{{\rm{d}}}^{3}{\boldsymbol{r}}^{\prime} ,\end{eqnarray} \tag{ E.5 }$$

with $r:= | {\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} |$ the distance from the point source to the point of interest. Working in prolate spheroidal coordinates we then have

$$\begin{eqnarray}&&\phi =-G\int \displaystyle \frac{1}{r}\rho (\xi ^{\prime} ,\eta ^{\prime} ){l}^{3}[\xi {{\prime} }^{2}-\eta {{\prime} }^{2}]{\rm{d}}\xi ^{\prime} {\rm{d}}\eta ^{\prime} {\rm{d}}\psi ^{\prime} .\end{eqnarray} \tag{ E.6 }$$

Following [114, 115], in prolate spheroidal coordinates, the ratio of l to r can be expressed as:

$$\begin{eqnarray}&&\displaystyle \frac{l}{r}=\displaystyle \sum _{n=0}^{\infty }(2n+1){P}_{n}(\eta ){P}_{n}(\eta ^{\prime} ){Q}_{n}(\xi ){P}_{n}(\xi ^{\prime} )\end{eqnarray} \tag{ E.7 }$$

$$\begin{eqnarray}&&+\ 2\displaystyle \sum _{n=1}^{\infty }(2n+1)\displaystyle \sum _{m=1}^{n}{\left(-1\right)}^{m}{\left(\displaystyle \frac{(n-m)!}{(n+m)!}\right)}^{2}\end{eqnarray} \tag{ E.8 }$$

$$\begin{eqnarray}&&\times {P}_{n}^{m}(\eta ){P}_{n}^{m}(\eta ^{\prime} ){Q}_{n}^{m}(\xi ){P}_{n}^{m}(\xi ^{\prime} )\cos \left[m(\psi -\psi ^{\prime} )\right],\end{eqnarray} \tag{ E.9 }$$

for $\xi \gt \xi ^{\prime}$, and

$$\begin{eqnarray}&&\displaystyle \frac{l}{r}=\displaystyle \sum _{n=0}^{\infty }(2n+1){P}_{n}(\eta ){P}_{n}(\eta ^{\prime} ){P}_{n}(\xi ){Q}_{n}(\xi ^{\prime} )+f\left(\cos \left[m\left(\psi -\psi ^{\prime} \right)\right]\right),\end{eqnarray} \tag{ E.10 }$$

for $\xi \lt \xi ^{\prime}$, where ${P}_{n}^{m}(x)$ and ${Q}_{n}^{m}(x)$ are the associated Legendre polynomials of the first and second kind and f is an unimportant function of $\cos [m(\psi -\psi ^{\prime} )]$ since, when inserting (E.10) into (E.6), this term, and the second term in (E.10), vanish once we integrate over ψ [114, 115]. For the prolate spheroid, we then end up with:

$$\begin{eqnarray}&&{\phi }_{{\rm{out}}}=-2\pi {l}^{2}G\displaystyle \sum _{n=0}^{\infty }(2n+1){Q}_{n}(\xi ){P}_{n}(\eta )\end{eqnarray} \tag{ E.11 }$$

$$\begin{eqnarray}&&\times {\int }_{-1}^{1}{\int }_{1}^{{\xi }_{0}}\rho (\xi ^{\prime} ,\eta ^{\prime} ){P}_{n}(\xi ^{\prime} ){P}_{n}(\eta ^{\prime} )\left[\xi {{\prime} }^{2}-\eta {{\prime} }^{2}\right]{\rm{d}}\xi ^{\prime} {\rm{d}}\eta ^{\prime} ,\end{eqnarray} \tag{ E.12 }$$

$$\begin{eqnarray}&&{\phi }_{{\rm{in}}}=-2\pi {l}^{2}G\displaystyle \sum _{n=0}^{\infty }(2n+1)\left[{Q}_{n}(\xi ){P}_{n}(\eta ){\int }_{-1}^{1}{\int }_{1}^{\xi }\rho (\xi ^{\prime} ,\eta ^{\prime} ){P}_{n}(\xi ^{\prime} ){P}_{n}(\eta ^{\prime} )\right.\end{eqnarray} \tag{ E.13 }$$

$$\begin{eqnarray}&&\left.+\ {P}_{n}(\xi ){P}_{n}(\eta ){\int }_{-1}^{1}{\int }_{\xi }^{{\xi }_{0}}\rho (\xi ^{\prime} ,\eta ^{\prime} ){Q}_{n}(\xi ^{\prime} ){P}_{n}(\eta ^{\prime} )\right][\xi {{\prime} }^{2}-\eta {{\prime} }^{2}]{\rm{d}}\xi ^{\prime} {\rm{d}}\eta ^{\prime} .\end{eqnarray} \tag{ E.14 }$$

Using (E.4) for $\rho (\xi ^{\prime} ,\eta ^{\prime} )$, we then find that:

$$\begin{eqnarray}&&{\phi }_{{\rm{out}}}^{{\rm{prolate}}}=-\displaystyle \frac{{GM}}{112l}\left(16{Q}_{0}({\xi }_{0}){F}_{1}(\xi ,\eta )+\displaystyle \frac{{\xi }^{2}-{\xi }_{0}^{2}}{{\xi }_{0}^{3}{\left[{\xi }_{0}^{2}-1\right]}^{2}}{F}_{2}(\xi ,\eta )\right.\end{eqnarray} \tag{ E.15 }$$

$$\begin{eqnarray}&&\left.+\ \displaystyle \frac{{\xi }^{2}({\xi }^{2}-1)}{{\xi }_{0}^{3}{\left[{\xi }_{0}^{2}-1\right]}^{2}}{F}_{3}(\xi ,\eta )\right),\end{eqnarray} \tag{ E.16 }$$

$$\begin{eqnarray}&&{\phi }_{{\rm{out}}}^{{\rm{prolate}}}=-\displaystyle \frac{{GM}}{7l}\left(7{Q}_{0}(\xi )-10{Q}_{2}(\xi ){P}_{2}(\eta )+3{Q}_{4}(\xi ){P}_{4}(\eta )\right),\end{eqnarray} \tag{ E.17 }$$

where

$$\begin{eqnarray}&&{F}_{1}(\xi ,\eta ):= 7-10{P}_{2}(\xi ){P}_{2}(\eta )+3{P}_{4}(\xi ){P}_{4}(\eta ),\end{eqnarray} \tag{ E.18 }$$

$$\begin{eqnarray}&&{F}_{2}(\xi ,\eta ):= 14\left[5{\xi }_{0}^{2}-7{\xi }_{0}^{4}+{\xi }^{2}(3{\xi }_{0}^{2}-1)\right]\end{eqnarray} \tag{ E.19 }$$

$$\begin{eqnarray}&&+\ 40{P}_{2}(\xi ){P}_{2}(\eta )\left[2{\xi }_{0}^{2}-4{\xi }_{0}^{4}+{\xi }^{2}(3{\xi }_{0}^{2}-4)+3{\xi }^{4}\right]\end{eqnarray} \tag{ E.20 }$$

$$\begin{eqnarray}&&-\ 8{P}_{4}(\xi ){P}_{4}(\eta )\left[10{\xi }_{0}^{2}-6{\xi }_{0}^{4}-{\xi }^{2}(13+6{\xi }_{0}^{2})+15{\xi }^{4}\right],\end{eqnarray} \tag{ E.21 }$$

$$\begin{eqnarray}&&{F}_{4}(\xi ,\eta ):= 60{P}_{2}(\eta )[{\xi }_{0}^{4}-6{\xi }_{0}^{2}+2{\xi }^{2}-3{\xi }^{4}]\end{eqnarray} \tag{ E.22 }$$

$$\begin{eqnarray}&&+\ 5{P}_{4}(\eta )(21{\xi }^{2}-11)[3{\xi }_{0}^{2}-{\xi }^{2}(1+7{\xi }_{0}^{2})+5{\xi }^{4}],\end{eqnarray} \tag{ E.23 }$$

$$\begin{eqnarray}&&{P}_{4}(x):= \displaystyle \frac{1}{8}(3-30{x}^{2}+35{x}^{4}),\end{eqnarray} \tag{ E.24 }$$

$$\begin{eqnarray}&&{Q}_{4}(x):= \displaystyle \frac{1}{48}\xi (110-210{\xi }^{2})+{P}_{4}(\xi ){Q}_{0}(\xi ),\end{eqnarray} \tag{ E.25 }$$

and we have used the orthogonality relationship of Legendre polynomials of the first kind:

$$\begin{eqnarray}&&{\int }_{-1}^{1}{P}_{n}(x){P}_{m}(x){\rm{d}}{x}=\displaystyle \frac{2}{2n+1}{\delta }_{{mn}}.\end{eqnarray} \tag{ E.26 }$$

In contrast to the uniform case, we now have Legendre polynomials of the fourth degree. Also note that, unlike in the spherical case, the potential outside the BEC spheroid is different to the uniform spheroid. To obtain the oblate potentials in oblate coordinates, we just make the changes $\xi \to {\rm{i}}\xi$, ${\xi }_{0}\to {\rm{i}}{\xi }_{0}$ and $l\to -{\rm{i}}{l}$.

It is also possible to find the potentials in cylindrical coordinates by taking the inverse transformations:

$$\begin{eqnarray}&&\xi =\displaystyle \frac{1}{2l}\left[\sqrt{{r}^{2}+{\left(z+l\right)}^{2}}+\sqrt{{r}^{2}+{\left(z-l\right)}^{2}}\right],\end{eqnarray} \tag{ E.27 }$$

$$\begin{eqnarray}&&\eta =\displaystyle \frac{1}{2l}\left[\sqrt{{r}^{2}+{\left(z+l\right)}^{2}}-\sqrt{{r}^{2}+{\left(z-l\right)}^{2}}\right]\end{eqnarray} \tag{ E.28 }$$

for prolate spheroidal coordinates, and (taking $l\to -{\rm{i}}{l}$):

$$\begin{eqnarray}&&\xi =\displaystyle \frac{1}{2l}\left[\sqrt{{r}^{2}+{\left(z-{\rm{i}}{l}\right)}^{2}}+\sqrt{{r}^{2}+{\left(z+{\rm{i}}{l}\right)}^{2}}\right]\equiv \displaystyle \frac{1}{\sqrt{2}l}{E}_{o},\end{eqnarray} \tag{ E.29 }$$

$$\begin{eqnarray}&&\eta =\displaystyle \frac{{\rm{i}}}{2l}\left[\sqrt{{r}^{2}+{\left(z-{\rm{i}}{l}\right)}^{2}}-\sqrt{{r}^{2}+{\left(z+{\rm{i}}{l}\right)}^{2}}\right]\equiv \displaystyle \frac{1}{\sqrt{2}l}{F}_{o},\end{eqnarray} \tag{ E.30 }$$

for oblate spheroidal coordinates, where E_o is defined in the previous section, and

$$\begin{eqnarray}&&{F}_{o}:= \sqrt{{l}^{2}-{r}^{2}-{z}^{2}+\sqrt{{\left({l}^{2}-{r}^{2}\right)}^{2}+2\left({l}^{2}+{r}^{2}\right){z}^{2}+{z}^{4}}}.\end{eqnarray} \tag{ E.31 }$$

In the appropriate limit, i.e. $a\to c\to R$ (and so $l\to 0$), it is possible to show that these potentials become the spherical BEC potentials provided in appendix B. To calculate E_G for the different spheroidal configurations, the procedure in the previous section can be followed with the uniform potentials and density functions replaced with those above. Alternatively, it is possible to integrate over spheroidal surfaces of constant density using the area element ${l}^{2}{\left[\left({\xi }^{2}-\eta {\left(\xi \right)}^{2}\right)\left(1-\eta {\left(\xi \right)}^{2}+\left({\xi }^{2}-1\right)\eta ^{\prime} {\left(\xi \right)}^{2}\right)\right]}^{1/2}$ for prolate coordinates where $\eta ^{\prime} := {\rm{d}}\eta (\xi )/{\rm{d}}\xi$.

Here we calculate the gravitational potential due to a spheroidal BEC in the Gaussian limit for small ellipticity values. We work in spherical coordinates to easily compare to the spherical BEC case. In general, the potential can be calculated from (E.6), which, in spherical coordinates is

$$\begin{eqnarray}&&\phi (r,\theta ,\psi )=-G{\int }_{0}^{2\pi }{\int }_{0}^{\pi }{\int }_{0}^{\infty }\displaystyle \frac{\rho (r^{\prime} ,\theta ^{\prime} ,\psi ^{\prime} )}{| {\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} | }r{{\prime} }^{2}\sin \theta ^{\prime} {\rm{d}}{r}^{\prime} {\rm{d}}\theta ^{\prime} {\rm{d}}\psi ^{\prime} .\end{eqnarray} \tag{ F.1 }$$

Since the spheroidal density does not depend on ψ, we can set ψ = 0 such that

$$\begin{eqnarray}&&| {\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} | =\sqrt{{r}^{2}+r{{\prime} }^{2}-2{rr}^{\prime} (\sin \theta \sin \theta ^{\prime} \cos \psi ^{\prime} +\cos \theta \cos \theta ^{\prime} )}.\end{eqnarray} \tag{ F.2 }$$

We can then expand the above in terms of Legendre polynomials of the first kind:

$$\begin{eqnarray}\displaystyle \frac{1}{2\pi }{\int }_{0}^{2\pi }\displaystyle \frac{1}{| {\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} | }{\rm{d}}\psi ^{\prime} =\left\{\begin{array}{cl}\displaystyle \frac{1}{r}\displaystyle \sum _{n=0,\infty }{\left(\displaystyle \frac{r^{\prime} }{r}\right)}^{n}{P}_{n}(\cos \theta ){P}_{n}(\cos \theta ^{\prime} ) & \mathrm{if}\ r\gt r^{\prime} ,\\ \displaystyle \frac{1}{r^{\prime} }\displaystyle \sum _{n=0,\infty }{\left(\displaystyle \frac{r}{r^{\prime} }\right)}^{n}{P}_{n}(\cos \theta ){P}_{n}(\cos \theta ^{\prime} ) & \mathrm{if}\ r\lt r^{\prime} ,\end{array}\right.\end{eqnarray} \tag{ F.3 }$$

so that (E.6) becomes

$$\begin{eqnarray}&&\phi (r,\theta )=-2\pi {{GP}}_{n}(\cos \theta ){\int }_{0}^{\pi }\left[\displaystyle \frac{1}{{r}^{n+1}}{\int }_{0}^{r}r{{\prime} }^{n+2}\right.\end{eqnarray} \tag{ F.4 }$$

$$\begin{eqnarray}&&\left.+\ {r}^{n}{\int }_{r}^{\infty }r{{\prime} }^{1-n}\right]\rho (r^{\prime} ,\theta ^{\prime} ){P}_{n}(\cos \theta ^{\prime} )\sin \theta ^{\prime} {\rm{d}}{r}^{\prime} {\rm{d}}\theta ^{\prime} .\end{eqnarray} \tag{ F.5 }$$

The density of the spheroidal BEC in the Gaussian limit is given by (43). Assuming the spheroid to have the same volume as a sphere with radius R, then in the limit of small ellipticity, the density function becomes

$$\begin{eqnarray}&&\rho (r,\theta )={\rho }_{{\rm{sphere}}}^{{\rm{Gauss}}}\left[1\pm \displaystyle \frac{2{r}^{2}{e}^{2}}{3{R}^{2}}{P}_{2}(\cos \theta )\right.\end{eqnarray} \tag{ F.6 }$$

$$\begin{eqnarray}&&-\displaystyle \frac{{r}^{2}{e}^{4}}{135{R}^{2}}\left(35-\displaystyle \frac{14{r}^{2}}{{R}^{2}}-10(\pm 7\beta +\displaystyle \frac{2{r}^{2}}{{R}^{2}}){P}_{2}(\cos \theta )\right.\end{eqnarray} \tag{ F.7 }$$

$$\begin{eqnarray}&&\left.\left.-\displaystyle \frac{36{r}^{2}}{{R}^{2}}{P}_{4}(\cos \theta )\right)\right],\end{eqnarray} \tag{ F.8 }$$

where ${\rho }_{{\rm{sphere}}}^{{\rm{Gauss}}}$ is given by (C.1); + is for the prolate case, which has β = 1; and − is for the oblate case, which has β = 2. Inserting the density expression into (F.4), and using the orthogonality relationship (E.26), we find

$$\begin{eqnarray}&&\phi (r,\theta )={\phi }_{{\rm{sphere}}}^{{\rm{Gauss}}}\pm \displaystyle \frac{{{GMe}}^{2}}{6{{Rr}}^{3}}{P}_{2}(\cos \theta )\left[e(r){X}_{2}(r)-21{R}^{5}\mathrm{erf}\left(\displaystyle \frac{r}{R}\right)\right]\end{eqnarray} \tag{ F.9 }$$

$$\begin{eqnarray}&&+\displaystyle \frac{{{GMe}}^{4}}{{R}^{3}}\left(\displaystyle \frac{1}{45}e(r)({r}^{2}+{R}^{2})\right.\end{eqnarray} \tag{ F.10 }$$

$$\begin{eqnarray}&&\pm \displaystyle \frac{1}{63\beta {r}^{3}}\left[e(r){X}_{4}(r)-21{R}^{5}\mathrm{erf}\left(\displaystyle \frac{r}{R}\right)\right]{P}_{2}(\cos \theta )\end{eqnarray} \tag{ F.11 }$$

$$\begin{eqnarray}&&\left.+\displaystyle \frac{1}{140{r}^{5}}\left[e(r){Y}_{4}(r)-105{R}^{7}\mathrm{erf}\left(\displaystyle \frac{r}{R}\right)\right]{P}_{4}(\cos \theta )\right),\end{eqnarray} \tag{ F.12 }$$

where

$$\begin{eqnarray}&&e(r):= \displaystyle \frac{2}{\sqrt{\pi }}{{\rm{e}}}^{-{r}^{2}/{R}^{2}},\end{eqnarray} \tag{ F.13 }$$

$$\begin{eqnarray}&&{X}_{2}(r):= 2{r}^{3}+3{{rR}}^{2},\end{eqnarray} \tag{ F.14 }$$

$$\begin{eqnarray}&&{X}_{4}(r):= \pm 2\beta {r}^{5}+14{r}^{3}{R}^{2}+21{{rR}}^{4},\end{eqnarray} \tag{ F.15 }$$

$$\begin{eqnarray}&&{Y}_{4}(r):= r(8{r}^{6}+28{r}^{4}{R}^{2}+70{r}^{2}{R}^{4}+105{R}^{6}),\end{eqnarray} \tag{ F.16 }$$

${\phi }_{{\rm{sphere}}}^{{\rm{Gauss}}}$ is given by (C.3); and, again, + is for the prolate case, which has β = 1; and − is for the oblate case, which has β = 2.