Detection of gravitational waves using parametric resonance in Bose–Einstein condensates

The recent detection of gravitational waves is revolutionizing our understanding of the Universe. Though our current observations have only detected black holes and neutron stars with LIGO/Virgo [1], it is expected that future gravitational wave detectors may also see phenomena such as supernovae, extreme-mass inspirals, and a stochastic gravitational wave background [2–4].

One of the challenges of gravitational wave astronomy is that current detectors are sensitive to only a small range of frequencies. Constructing a multitude of detectors (both terrestrial and in space) is thus required to cover the entire gravitational wave spectrum [5, 6]. Much effort has been exerted into designing detectors for sub-kHz sources, such as the laser interferometer space antenna (0.1 mHz–1 Hz) [7], the Einstein telescope (1–10 000 Hz) [8, 9] and pulsar timing arrays like the International pulsar timing array and European pulsar timing array (∼nHz) [10, 11]. Atom interferometers have also been proposed to detect gravitational waves in the millihertz to decahertz regimes [12–14]. Though the Einstein telescope is proposed to be sensitive up to 10 kHz, its optimal sensitivity will occur around the same frequency range as LIGO. Additional gravitational wave detector technologies have also been proposed, including a superfluid detector for continuous gravitational waves (0.1–1.5 kHz) [15] and a satellite-based cold atom interferometer for mid-frequency gravitational waves (30 mHz–10 Hz) [16].

In the frequency regime probed by LIGO, Virgo, and KAGRA (above 100 Hz and below 1 kHz), the most common sources are binary black holes, binary neutron stars, and black hole–neutron star mergers [17]. In the kHz regime, it is predicted that transient sources include lower-mass black holes and neutron star mergers at frequencies outside the range that LIGO can currently observe (up to several kHz) [18–20] and magnetars (0.5–2 kHz) [21]. Depending on the model, there could potentially even be primordial black holes in the kHz domain [22]. In the continuous regime, there may be neutron stars/pulsars (tens-hundreds of Hz) [23, 24] and boson clouds (extending into the kHz regime and above) [25].

One intriguing proposal [26] suggests using a uniform (non-rotating) Bose–Einstein condensate (BEC) as a gravitational wave detector. By noting that an incoming gravitational wave is able to create phonons within the condensate and comparing the state of the phonons before and after the gravitational wave interacts with the BEC, it may be possible to gain knowledge about the amplitude of the gravitational wave. To estimate the classical amplitude ⁵ the quantum Fisher information H , which describes how much quantum information can be obtained from a single measurement of an arbitrary parameter in a quantum system. By doing M measurements of the amplitude of the gravitational wave (for example, using quantum dots or phonon evaporation [26]), the quantum Cramer–Rao bound states that the resultant uncertainty in the measurement of is $\langle {({\Delta}{\epsilon})}^{2}\rangle \geqslant \frac{1}{M{H}_{{\epsilon}}}$ [29]. Increasing either the number of measurements of the system or the quantum Fisher information obtained from each measurement therefore improves sensitivity to gravitational wave detection.

Such strategies of using cold atoms and BECs to investigate fundamental physics and curved spacetime have increased in prominence in recent years. For instance, it has been shown that BECs can be used in the study of analogue gravity [30, 31] to simulate different aspects of curved spacetime, such as black holes [32], inflation [33, 34], false vacuum decay [35], Hawking radiation [36, 37], as well as gravitational waves [38, 39]. In addition to gravitational wave detection [12–14, 16], cold atom interferometers are also expected to provide new insights and measurements into the gravitational constant [40], the equivalence principle [41] and dark energy [42, 43].

However, this proposal of using a BEC to detect gravitational waves faces significant challenges. It has been argued [44] that using a BEC to detect gravitational waves would be quite difficult as LIGO is able to use large numbers (e.g. ratio of interferometric arm length to laser wavelength, photon number in the interferometer), to achieve the requisite sensitivity to gravitational waves. By contrast, a BEC is several orders of magnitude smaller in size and its number of particles is many orders of magnitude smaller than the number of photons in LIGO. In addition, it may be better to use a non-uniform condensate to detect gravitational waves [44, 45]. In response it has been has argued that such large numbers may be achieved partly through the number of atoms in the condensate and the number of squeezed phonons present [46]. We emphasize that initially squeezing is key to the original proposal of [26]. In [47], we also pointed out that, while a BEC could in principle be used to detect gravitational waves, it is necessary to squeeze the initial state of the BEC phonons well beyond current experimental capabilities. Furthermore, the natural decay of the BEC by Beliaev damping limits the sensitivity, especially at the higher frequencies of interest. Decoherence and the decay of sensitivity were considered in more detail in [48].

We investigate here the possibility of improving the sensitivity of a BEC detector by modulating the speed of sound to induce a parametric resonance with the BEC phonons, which then boosts their sensitivity to gravitational waves. In contrast to previous studies [26, 47], we find that parametric resonance causes the BEC to be more sensitive to gravitational waves at lower frequencies than at higher frequencies, though detection of higher frequency gravitational waves may still be possible. The optimum sensitivity depends upon the parameters of the condensate such as its length, speed of sound and phonon frequency.

Noting that LIGO detects gravitational waves by exploiting ratios of large numbers in order to improve the sensitivity by 20 orders of magnitude [44], in our proposal close to 20 orders of magnitude can be obtained through the parametric amplification of phonons by modulations of the speed of sound. These modulations induce instabilities in the phononic modes, which we then exploit. In addition, by repeatedly measuring the condensate over the course of a year, we exploit the ratio between the duration of a year and the period of the gravitational wave. Approximately 3–4 orders of magnitude in sensitivity can be obtained by running the detector over the course of a year. The remaining sensitivity is obtained through the parametric resonance effects (with the amount of sensitivity determined by the experimental parameters of the condensate, such as the mass of the atoms, lifetime of the BEC, size of the BEC, and the speed of sound).

It should be emphasized that, in contrast to [26, 47], the initial squeezing of phonons is not a requirement of our proposal. Instead, gravitational wave detection could potentially be achieved through increased laser power, larger condensates, using a similar suspension system as LIGO to minimize vibrations in the optical setup, and squeezed light generating the trap.

We emphasize that we will be modulating the speed of sound, rather than the trap parameters, as trap parameter modulations would be a direct source of phonons. For simplicity, we will be assuming a cubic condensate of uniform density, though the same analysis can be applied for more realistic trap geometries. They key to the proposal is to choose a phonon mode on resonance with the gravitational wave. If such a mode exists, then a BEC could potentially be used to the detect a gravitational wave, regardless of trap geometry. The trap geometry only affects the sensitivity to that wave (though any loss in sensitivity can be compensated by adjusting other parameters of the trap).

The outline of the subsequent sections is as follows: in section 2, we review the theory behind BECs in curved spacetime and solve the equations of motion of the phonons as well as the Bogoliubov coefficients. We also discuss the effect of parametric resonance and non-linearities in our model. In section 3, we then show how quantum Fisher information is used to determine the sensitivity to gravitational waves. Next, in section 4, we discuss how Beliaev damping will limit the sensitivity of the BEC to gravitational waves at high frequencies. In section 5, we discuss the possible sources that a BEC could observe as well as consider the sensitivity of a BEC experiment to gravitational waves, comment on expected noise sources, and discuss what may be necessary to make this proposal a reality. Section 6 summarizes our conclusions and presents possible directions for future work.

2.1. Equation of motion

We first derive the equations of motion and Bogoliubov coefficients for BEC phonons in a curved spacetime. Note that similar derivations have also been done in [26, 31, 49–51]. The Lagrangian for a BEC in a curved spacetime is

$$\begin{equation}\mathcal{L}={g}^{\mu \nu }{\partial }_{\mu }\phi {\partial }_{\nu }{\phi }^{{\dagger}}-{m}^{2}\vert \phi {\vert }^{2}-U(\vert \phi {\vert }^{2}),\end{equation} \tag{ 1 }$$

where g_μν is the metric, ϕ is the scalar field of the bosonic atoms, m is the mass of the atoms, and U(|ϕ|²) = λ(t)|ϕ|⁴ > 0 is the interaction potential, and λ(t) is the interaction strength (which we relabel as λ for the duration of this paper).

Given the definition $\phi =\hat{\phi }{\text{e}}^{\text{i}\chi }$ with real $\hat{\phi }$ and χ, we will assume a slowly varying amplitude, i.e. that ${(\partial \hat{\phi })}^{2}\ll {\hat{\phi }}^{2}{(\partial \chi )}^{2}$. This corresponds to the non-relativistic approximation (as discussed in the appendix of [47]) of the BEC, as the phase of ${\text{e}}^{-\text{i}\frac{m{c}^{2}t}{\hslash }}$ has not yet been separated from the scalar field in the Lagrangian. We can extremize $\mathcal{L}$ with respect to $\hat{\phi }$ to find:

$$\begin{equation}{\hat{\phi }}^{2}=\frac{1}{2\lambda }\left[{\partial }^{\mu }\chi {\partial }_{\mu }\chi -{m}^{2}\right],\end{equation} \tag{ 2 }$$

and so

$$\begin{equation}\mathcal{L}=\frac{1}{2\lambda }{\left[{g}^{\mu \nu }{\partial }_{\mu }\chi {\partial }_{\nu }\chi -{m}^{2}\right]}^{2}.\end{equation} \tag{ 3 }$$

We note that at the boundaries of the trap, there will not be a slowly varying-amplitude, but rather a sharp cutoff. However, we will neglect this effect as we are assuming an ideal system in which the density of the trap is constant throughout, except right at the boundaries. The phonons cannot propagate outside the walls, and thus we can approximately use Dirichlet boundary conditions. It should be noted that Bose–Einstein condensation can occur when using other boundary conditions (see, for example, [52]), though we focus on Dirichlet boundary conditions as an idealization.

Let $\chi =\frac{m{c}^{2}t}{\hslash }+f(t)+\pi ({x}^{\mu })$, where $\pi ({x}^{\mu })\in \mathfrak{R}$ is a Goldstone boson representing the acoustic perturbation (phonon) field of the condensate. These phonons can be represented as fluctuations in the phase [53]. Inserting this ansatz into equation (3), we get

$$\begin{equation}S=\int \frac{{\mathrm{d}}^{4}x}{4\lambda }\sqrt{-g}{\left\{(\dot{f}{\delta }_{0}^{\nu }+{\partial }^{\nu }\pi )(\dot{f}{\delta }_{0}^{\mu }+{\partial }^{\mu }\pi ){g}_{\mu \nu }-{m}^{2}\right\}}^{2}.\end{equation} \tag{ 4 }$$

Let us take g_μν = η_μν + h_μν , where [h_μν ] ≪ 1 describes the distortion of spacetime due to a gravitational wave. As we showed in [47], the action can be written as

$$\begin{equation}S\approx \int \frac{{\mathrm{d}}^{4}x}{4\lambda }\left[\vert \dot{\pi }{\vert }^{2}(6{\dot{f}}^{2}-2{m}^{2})+(2{\dot{f}}^{2}-2{m}^{2})({\eta }_{ij}+{h}_{ij}){\partial }^{i}\pi {\partial }^{j}\pi \right],\end{equation} \tag{ 5 }$$

where we are working in the (+, −, −, −) convention and have expanded to second order in π, assuming that the higher-order terms in π can be neglected and the first-order term integrates to zero on the boundary. Calculating the Euler–Lagrange equation gives

$$\begin{equation}\ddot{\pi }+\frac{4\lambda \dot{\pi }}{6{\dot{f}}^{2}-2{m}^{2}}\frac{\partial }{\partial t}\left[\frac{6{\dot{f}}^{2}-2{m}^{2}}{4\lambda }\right]+\frac{{\dot{f}}^{2}-{m}^{2}}{3{\dot{f}}^{2}-{m}^{2}}\left({\eta }_{ij}+{h}_{ij}\right){\partial }^{i}{\partial }^{j}\pi =0.\end{equation} \tag{ 6 }$$

We will now consider a BEC with its scattering length modulated using a Feshbach resonance, such that it will have an oscillating speed of sound. Let $\frac{{\dot{f}}^{2}-{m}^{2}}{3{\dot{f}}^{2}-{m}^{2}}\equiv {c}_{s}^{2}(t)$, where f(t) is a function that depends on the interaction potential λ. In this paper, we will consider modulations of the speed of sound of the form ${c}_{s}^{2}(t)={\bar{c}}_{s}^{2}\left(1+2a\,\mathrm{sin}\,{{\Omega}}_{B}t\right)$, where ${\bar{c}}_{s}$ is the average speed of sound, Ω_B is the oscillation frequency of an applied magnetic field that changes the speed of sound by changing λ, and a ≪ 1. Assuming for simplicity that $\frac{4\lambda \dot{\pi }}{6{\dot{f}}^{2}-2{m}^{2}}\frac{\partial }{\partial t}\left[\frac{6{\dot{f}}^{2}-2{m}^{2}}{4\lambda }\right]$ is negligible ⁶ , we find the equation of motion

$$\begin{equation}\ddot{\pi }+{\bar{c}}_{s}^{2}(1+2a\,\mathrm{sin}\,{{\Omega}}_{B}t)({\eta }_{ij}+{h}_{ij}){\partial }^{i}{\partial }^{j}\pi =0.\end{equation} \tag{ 7 }$$

2.2. Bogoliubov coefficients

Let $\pi ={\text{e}}^{\text{i}\vec{k}\cdot \bar{x}}\psi (t)$ and assume a linear dispersion relation. We will also work in the TT-gauge and assume no cross-polarization for simplicity,

$$\begin{equation}{h}_{\mu \nu }=\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {h}_{+}(t)\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -{h}_{+}(t)\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right).\end{equation} \tag{ 8 }$$

Taking h₊(t) =  sin Ωt, where is the amplitude of the continuous gravitational wave and Ω is its frequency, the equation of motion (7) becomes

$$\begin{equation}\ddot{\psi }+{\omega }^{2}(1+2a\,\mathrm{sin}\,{{\Omega}}_{B}t)(1+\tilde{{\epsilon}}\,\mathrm{sin}\,{\Omega}t)\psi =0,\end{equation} \tag{ 9 }$$

where ${\omega }^{2}={\bar{c}}_{s}^{2}{k}^{2}$ and $\tilde{{\epsilon}}=\frac{({k}_{x}^{2}-{k}_{y}^{2})}{\vert k{\vert }^{2}}{\epsilon}$, which we see can be interpreted as having a modulating speed of sound $\ddot{\psi }+{k}^{2}{c}_{s}{(t)}^{2}\psi =0$. Part of the modulation arises from Feshbach resonance, with the other part of the modulation arising due to gravitational waves.

Let ${\omega }_{1}{(t)}^{2}={\omega }^{2}(1+\tilde{{\epsilon}}\,\mathrm{sin}\,{\Omega}t+2a\,\mathrm{sin}\,{{\Omega}}_{B}t+2a\tilde{{\epsilon}}\,\mathrm{sin}\,{\Omega}t\,\mathrm{sin}\,{{\Omega}}_{B}t)$. As we want to consider the gravitational wave on resonance with the frequency of the modulations of the speed of sound, we will take Ω_B = Ω. We can expand to first-order in and second-order in a (neglecting the $\mathcal{O}(a\tilde{{\epsilon}})$ terms as $a\tilde{{\epsilon}}\ll a$). We will see that it is the presence of the $\mathcal{O}(a)$ terms that induce parametric resonance within the condensate.

To solve equation (9) in this approximation, let us assume that we can decompose the solutions as [54]

$$\begin{equation}\pi =\alpha (t)\frac{{\text{e}}^{-\text{i}\int {\omega }_{1}({t}^{\prime })\mathrm{d}{t}^{\prime }}}{\sqrt{2{\omega }_{1}(t)}}+\beta (t)\frac{{\text{e}}^{\text{i}\int {\omega }_{1}({t}^{\prime })\mathrm{d}{t}^{\prime }}}{\sqrt{2{\omega }_{1}(t)}},\end{equation} \tag{ 10 }$$

$$\begin{equation}\dot{\pi }=-i\alpha (t)\sqrt{\frac{{\omega }_{1}(t)}{2}}{\text{e}}^{-\text{i}\int {\omega }_{1}({t}^{\prime })\mathrm{d}{t}^{\prime }}+i\beta (t)\sqrt{\frac{{\omega }_{1}(t)}{2}}{\text{e}}^{\text{i}\int {\omega }_{1}({t}^{\prime })\mathrm{d}{t}^{\prime }},\end{equation} \tag{ 11 }$$

where α(t) and β(t) are time-dependent Bogoliubov coefficients that satisfy the coupled differential equations

$$\begin{equation}\dot{\alpha }=\frac{{\dot{\omega }}_{1}}{2{\omega }_{1}}{e}^{2i\int {\omega }_{1}({t}^{\prime })\mathrm{d}{t}^{\prime }}\beta (t)\qquad \dot{\beta }=\frac{{\dot{\omega }}_{1}}{2{\omega }_{1}}{e}^{-2i\int {\omega }_{1}({t}^{\prime })\mathrm{d}{t}^{\prime }}\alpha (t).\end{equation} \tag{ 12 }$$

We will solve these coupled equations using the method outlined in [55]. Setting $\hat{\alpha }={\text{e}}^{-\text{i}\bar{\omega }t}\alpha$ and $\hat{\beta }={\text{e}}^{\text{i}\bar{\omega }t}\beta$, where $\bar{\omega }=\frac{1}{T}{\int }_{0}^{T}{\omega }_{1}(t)\mathrm{d}t$, we obtain from (12),

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{matrix}\hfill \hat{\alpha }\hfill \\ \hfill \hat{\beta }\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill -i\bar{\omega }\hfill & \hfill \frac{{\dot{\omega }}_{1}}{2{\omega }_{1}}{e}^{2i\int \delta \omega ({t}^{\prime })\mathrm{d}{t}^{\prime }}\hfill \\ \hfill \frac{{\dot{\omega }}_{1}}{2{\omega }_{1}}{e}^{-2i\int \delta \omega ({t}^{\prime })\mathrm{d}{t}^{\prime }}\hfill & \hfill i\bar{\omega }\hfill \end{matrix}\right)\left(\begin{matrix}\hfill \hat{\alpha }\hfill \\ \hfill \hat{\beta }\hfill \end{matrix}\right),\end{equation} \tag{ 13 }$$

where $\delta \omega (t)={\omega }_{1}(t)-\bar{\omega }$. Let $\hat{\alpha }={\hat{\alpha }}_{0}+b{\hat{\alpha }}_{1}+{b}^{2}{\hat{\alpha }}_{2}$ and $\hat{\beta }={\hat{\beta }}_{0}+b{\hat{\beta }}_{1}+{b}^{2}{\hat{\beta }}_{2}$, where b = 2a + ≪ 1. Note that, up to second order in b, $\bar{\omega }=\omega \left(1-\frac{{b}^{2}}{16}\right)$.

Substitution into equation (13), solving order-by-order, and restricting ourselves to solving the system after each period, we can define:

$$\begin{equation}\left(\begin{matrix}\hfill \hat{\alpha }(T)\hfill \\ \hfill \hat{\beta }(T)\hfill \end{matrix}\right)=M\left(\begin{matrix}\hfill \hat{\alpha }(0)\hfill \\ \hfill \hat{\beta }(0)\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill {M}_{11}\hfill & \hfill {M}_{12}\hfill \\ \hfill {M}_{12}^{\ast }\hfill & \hfill {M}_{22}\hfill \end{matrix}\right)\left(\begin{matrix}\hfill \hat{\alpha }(0)\hfill \\ \hfill \hat{\beta }(0)\hfill \end{matrix}\right),\end{equation} \tag{ 14 }$$

where after a period T = 2π/Ω, we have

$$\begin{equation}{M}_{11}=\frac{{e}^{-2i\pi q}\left(4+{q}^{2}\left({b}^{2}\left(-\left(-8i\pi {q}^{3}+2i\pi q+{e}^{4i\pi q}-1\right)\right)+64{q}^{2}-32\right)\right)}{4{\left(1-4{q}^{2}\right)}^{2}},\end{equation} \tag{ 15 }$$

$$\begin{equation}{M}_{12}=\frac{bq\left(-3ibq+4{q}^{2}-4\right)\mathrm{sin}(2\pi q)}{4\left(4{q}^{4}-5{q}^{2}+1\right)},\end{equation} \tag{ 16 }$$

$$\begin{equation}{M}_{22}=\frac{{e}^{-2i\pi q}\left(-{b}^{2}{q}^{2}+{e}^{4i\pi q}\left(4+{q}^{2}\left({b}^{2}\left(-8i\pi {q}^{3}+2i\pi q+1\right)+64{q}^{2}-32\right)\right)\right)}{4{\left(1-4{q}^{2}\right)}^{2}}.\end{equation} \tag{ 17 }$$

We have assumed that ${\hat{\alpha }}_{1}(0)={\hat{\beta }}_{1}(0)={\hat{\alpha }}_{2}(0)={\hat{\beta }}_{2}(0)=0$. We see that the 1st order resonance (i.e. to 2nd order in b) occurs when $q\approx \frac{1}{2}$ (see figure 1).

Zoom In Zoom Out Reset image size

After N periods of the gravitational wave oscillations (or equivalently after N periods of the speed of sound modulations), we have

$$\begin{equation}\left(\begin{matrix}\hfill \hat{\alpha }(NT)\hfill \\ \hfill \hat{\beta }(NT)\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill {e}^{-2\pi iqN}\alpha (NT)\hfill \\ \hfill {e}^{2\pi iqN}\beta (NT)\hfill \end{matrix}\right)={M}^{N}\left(\begin{matrix}\hfill \hat{\alpha }(0)\hfill \\ \hfill \hat{\beta }(0)\hfill \end{matrix}\right).\end{equation} \tag{ 18 }$$

Diagonalizing M and writing $\left(\begin{matrix}\hfill \hat{\alpha }(0)\hfill \\ \hfill \hat{\beta }(0)\hfill \end{matrix}\right)={k}_{1}{\lambda }_{1}{\vec{x}}_{1}+{k}_{2}{\lambda }_{2}{\vec{x}}_{2}$, where ${\lambda }_{i},{\vec{x}}_{i}$ are the eigenvalues/eigenvectors of M and k_i are coefficients related to the eigenvalues, we see

$$\begin{equation}\left(\begin{matrix}\hfill \hat{\alpha }(NT)\hfill \\ \hfill \hat{\beta }(NT)\hfill \end{matrix}\right)={k}_{1}{\lambda }_{1}^{N}{\vec{x}}_{1}+{k}_{2}{\lambda }_{2}^{N}{\vec{x}}_{2}.\end{equation} \tag{ 19 }$$

Let λ₂ correspond to the eigenvalue whose magnitude is less than one around resonance. We can neglect this term because it will become negligible after many periods. Therefore, the only relevant quantities are

$$\begin{equation}{k}_{1}=-\frac{b\left(128{\zeta }^{3}-64i{\zeta }^{2}+3\zeta +i\right)+128{\zeta }^{2}-2}{4{\left(1-64{\zeta }^{2}\right)}^{3/2}}\end{equation} \tag{ 20 }$$

$$\begin{equation}{\lambda }_{1}=\pi {b}^{2}\left(\frac{5\sqrt{1-64{\zeta }^{2}}\zeta }{512{\zeta }^{2}-8}+\pi \left(2{\zeta }^{2}-\frac{1}{32}\right)\right)-\frac{1}{4}\pi b\sqrt{1-64{\zeta }^{2}}-1\end{equation} \tag{ 21 }$$

where $\zeta \equiv (q-\frac{1}{2})/b$, and

$$\begin{equation}{\vec{x}}_{1}=\left(\begin{matrix}\hfill {x}_{0}+b{x}_{1}+{b}^{2}{x}_{2}\hfill \\ \hfill 1\hfill \end{matrix}\right)\end{equation} \tag{ 22 }$$

with

$$\begin{equation}{x}_{0}=\sqrt{1-64{\zeta }^{2}}-8i\zeta \end{equation} \tag{ 23 }$$

$$\begin{equation}{x}_{1}=8i{\zeta }^{2}+\frac{\zeta (3+64\zeta (2\zeta -i))+i}{2\sqrt{1-64{\zeta }^{2}}}+4\zeta +\frac{3i}{16}\end{equation} \tag{ 24 }$$

$$\begin{equation}{x}_{2}=-16i{\zeta }^{3}+\frac{28{\zeta }^{2}}{3}-\frac{1}{12}i{\pi }^{2}\zeta +\frac{37i\zeta }{16}-\frac{3}{32}\end{equation} \tag{ 25 }$$

$$\begin{equation}+\frac{-219+128\zeta \left(4\zeta \left(2{\pi }^{2}\left(64{\zeta }^{2}-1\right)+16\zeta (4\zeta (-63+32\zeta (12\zeta +7i))-43i)+75\right)+29i\right)}{1536{\left(1-64{\zeta }^{2}\right)}^{3/2}}.\end{equation} \tag{ 26 }$$

Let us further define ${\alpha }_{2a+\tilde{{\epsilon}}}={\alpha }_{2a}+\tilde{{\epsilon}}{\alpha }_{\tilde{{\epsilon}}}$ and ${\beta }_{2a+\tilde{{\epsilon}}}={\beta }_{2a}+\tilde{{\epsilon}}{\beta }_{\tilde{{\epsilon}}}$. The goal is to calculate the Bogoliubov coefficients ${\alpha }_{\tilde{{\epsilon}}}$ and ${\beta }_{\tilde{{\epsilon}}}$, which correspond to the effect solely due to the gravitational wave interacting with the speed of sound modulation. Using the definition of the Bogoliubov coefficients, we can write ${\hat{a}}_{0}={\alpha }_{2a+\tilde{{\epsilon}}}{\hat{a}}_{2a+\tilde{{\epsilon}}}+{\beta }_{2a+\tilde{{\epsilon}}}^{\ast }{\hat{a}}_{2a+\tilde{{\epsilon}}}^{{\dagger}}$, where ${\hat{a}}_{0}$ is the annihilation operation in the mode decomposition with no gravitational wave or modulation and ${\hat{a}}_{2a+\tilde{{\epsilon}}}$ is the annihilation operator in the mode decomposition containing both a gravitational wave and modulation. Thus ${\hat{a}}_{2a+\tilde{{\epsilon}}}={\alpha }_{2a+\tilde{{\epsilon}}}^{\ast }{\hat{a}}_{0}-{\beta }_{2a+\tilde{{\epsilon}}}^{\ast }{\hat{a}}_{0}^{{\dagger}}$. We can also write ${\hat{a}}_{2a+\tilde{{\epsilon}}}={\alpha }_{\text{rel}}{\hat{a}}_{2a}-{\beta }_{\text{rel}}^{\ast }{\hat{a}}_{2a}^{{\dagger}}$. Therefore,

$$\begin{equation}{\hat{a}}_{2a+\tilde{{\epsilon}}}={\alpha }_{2a+\tilde{{\epsilon}}}^{\ast }\left({\alpha }_{2a}{\hat{a}}_{2a}+{\beta }_{2a}^{\ast }{\hat{a}}_{2a}^{{\dagger}}\right)-{\beta }_{2a+\tilde{{\epsilon}}}^{\ast }\left({\alpha }_{2a}^{\ast }{\hat{a}}_{2a}^{{\dagger}}+{\beta }_{2a}{\hat{a}}_{2a}\right)\end{equation} \tag{ 27 }$$

$$\begin{equation}=\left({\alpha }_{2a+\tilde{{\epsilon}}}^{\ast }{\alpha }_{2a}-{\beta }_{2a+\tilde{{\epsilon}}}^{\ast }{\beta }_{2a}\right){\hat{a}}_{2a}+\left({\alpha }_{2a+\tilde{{\epsilon}}}^{\ast }{\beta }_{2a}^{\ast }-{\beta }_{2a+\tilde{{\epsilon}}}^{\ast }{\alpha }_{2a}^{\ast }\right){\hat{a}}_{2a}^{{\dagger}}.\end{equation} \tag{ 28 }$$

With ${\alpha }_{2a+\tilde{{\epsilon}}}={\alpha }_{2a}+\tilde{{\epsilon}}{\alpha }_{\tilde{{\epsilon}}}$ and ${\beta }_{2a+\tilde{{\epsilon}}}={\beta }_{2a}+\tilde{{\epsilon}}{\beta }_{\tilde{{\epsilon}}}$, we see

$$\begin{equation}\ {\alpha }_{\text{rel}}^{\ast }=1+\tilde{{\epsilon}}\left({\alpha }_{\tilde{{\epsilon}}}^{\ast }{\alpha }_{2a}-{\beta }_{\tilde{{\epsilon}}}^{\ast }{\beta }_{2a}\right),\end{equation} \tag{ 29 }$$

$$\begin{equation}{\beta }_{\text{rel}}^{\ast }=-\tilde{{\epsilon}}\left({\alpha }_{\tilde{{\epsilon}}}^{\ast }{\beta }_{2a}^{\ast }-{\beta }_{\tilde{{\epsilon}}}^{\ast }{\alpha }_{2a}^{\ast }\right).\end{equation} \tag{ 30 }$$

2.3. Parametric amplification

2.3.1. Resonance

Equation (9) is a perturbed form of the Mathieu differential equation (as $a\tilde{{\epsilon}}$ is a sub-dominant term), which implies that our system includes parametric resonance effects [56, 57]. In figure 1(a), we show the locations of the first four resonances of the system, where the shading illustrates the instability bands associated with each individual resonance in the Mathieu equation (figure 1(a)) [58]. For a system lying within these bands (indicated by the shaded regions), evolution is unstable and will become non-perturbative after a certain amount of time, with the amount of time before non-perturbative effects occurs depending on the precise location within these regions. In figure 1(b), we see that if $q=\frac{\omega }{{\Omega}}$ is too far off the resonance, the instability band is exited, which limits production of particles and nonlinearities.

By considering phonons with frequencies that are within the instability bands (locations of parametric resonance) for a time t, we can enhance the sensitivity of the BEC to gravitational waves. In section 2.3.2, we address the maximum time that we can lie within the regions of instability before the system becomes non-perturbative in a.

We have focussed on parametric resonance, though a BEC could also include additional resonance effects, such as those arising from direct driving and mode coupling [59]. Mode–mode coupling is a negligible effect as it does not couple the Bogoliubov coefficients of our phonons. In our scheme of modulating the speed of sound, direct driving of phonons is a negligible, second-order effect. We want to minimize the number of extra phonons in the experiment arising from non-gravitational wave effects.

We note that we have approximated the BEC as homogeneous in a cubic trap. At the boundaries, this approximation may break down (where the level of homogeneity depends on the 'sharpness' of the corners of the trap). In this case, direct driving effects may affect the sensitivity to gravitational waves, though as a first-approximation, we are neglecting this effect as we are assuming that the corners of our box trap are perfectly sharp.

2.3.2. Non-linearities

Similar to [47], we have the constraint of $\dot{\pi }\ll \mu =m{c}_{s}^{2}$. By squaring both sides (as $\langle \dot{\pi }\rangle =0$) and writing $\dot{\pi }$ as a sum of creation and annihilation operators, we find that this is equivalent to ${\sum }_{\vec{k}}\vert {\beta }_{2a,\vec{k}}{\vert }^{2}\hslash {\omega }_{\vec{k}}\ll nm{c}_{s}^{2}{L}^{3}$, where L³ is the volume of the condensate (taking the BEC to be cubic with sides of length L), n is the number density, m is the mass of the atoms, and $\vert {\beta }_{2a,\vec{k}}{\vert }^{2}$ is the number of phonons present in each mode $\vec{k}$. ⁷

Rewriting $\vert {\beta }_{2a,\vec{k}}{\vert }^{2}$ in terms of $q=\frac{\omega }{{\Omega}}$ and assuming that the BEC is constructed such that only a single phonon mode dominates the resonance instability (i.e., the most unstable mode), we note that at each value of q, there is a maximum number of trap oscillations N(q) satisfying $\frac{\vert {\beta }_{2a,q}{\vert }^{2}\hslash {\omega }_{q}}{\rho {L}^{3}}\lessapprox 0.05$, where we relabel $\vert {\beta }_{2a,\vec{k}}{\vert }^{2}$ as |β_2a,q |² to explicitly indicate the dependence on q. Physically, N(q) describes the maximum time that a gravitational wave could be observed by the condensate before the condensate becomes non-linear (because, by assumption, we take the frequency of the speed of sound oscillation to be the same as the period of the gravitational wave). Noting that |β_2a,q |² = |k₁|²|λ₁|^2N , we see that ⁸

$$\begin{equation}N\leqslant \frac{\mathrm{log}\left[\frac{0.05mn{\bar{c}}_{s}^{2}{L}^{3}}{\pi \hslash q{f}_{\text{GW}}\vert {k}_{1}{\vert }^{2}}\right]}{2\,\mathrm{log}[\vert {\lambda }_{1}\vert ]}.\end{equation} \tag{ 31 }$$

In figure 2, we demonstrate the dependence of N on q, the size of the condensate, and the frequency of the speed of sound modulation. We consider a cross-section of this plot in figure 3. Throughout the remainder of this paper, all calculations and discussions will assume that this non-linearity condition is saturated.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We see in the above equation that a greater number density will increase the maximum value of N. However, this equation was derived by assuming a constant density. If the density were non-uniform, the maximum amplification time would be different for different phonon modes. The key is to generate a condensate with the appropriate parameters to maximize the value of N for the phonon mode of interest.

As phonons are affected by the gravitational wave, the gravitational wave's amplitude is imprinted on the phonon's density matrix. The question then becomes how to extract the gravitational wave amplitude from the density matrix. This can be done by exploiting techniques in quantum metrology, which allows measurements to be done on quantities in quantum systems that are not operator observables [60–63]. In this section, we will be determine the sensitivity of an undamped condensate at zero-temperature with no extra sources of noise in the system. In section 4, we will consider the more realistic case by investigating the sensitivity of a damped BEC.

Given a parameter $\tilde{{\epsilon}}$ in a quantum system, the quantum Cramer–Rao bound is [29]

$$\begin{equation}\langle {({\Delta}\tilde{{\epsilon}})}^{2}\rangle \geqslant \frac{1}{M{H}_{\tilde{{\epsilon}}}},\end{equation} \tag{ 32 }$$

where $\langle {({\Delta}\tilde{{\epsilon}})}^{2}\rangle$ is the expectation value of the uncertainty in $\tilde{{\epsilon}}$, M is the number of measurements of the system, and ${H}_{\tilde{{\epsilon}}}$ is the quantum Fisher information [29]

$$\begin{equation}{H}_{\tilde{{\epsilon}}}=\frac{8\left(1-\sqrt{F({\rho }_{\tilde{{\epsilon}}},{\rho }_{\tilde{{\epsilon}}+d\tilde{{\epsilon}}})}\right)}{d{\tilde{{\epsilon}}}^{2}},\end{equation} \tag{ 33 }$$

with $F({\rho }_{\tilde{{\epsilon}}},{\rho }_{\tilde{{\epsilon}}+d\tilde{{\epsilon}}})={\left[\mathrm{T}\mathrm{r}\,\sqrt{{\rho }_{\tilde{{\epsilon}}}\sqrt{{\rho }_{\tilde{{\epsilon}}+d\tilde{{\epsilon}}}}{\rho }_{\tilde{{\epsilon}}}}\right]}^{2}$ quantifying the overlap (fidelity) between the states ${\rho }_{\tilde{{\epsilon}}}$ and ρ_+d [64, 65]. We find the sensitivity to gravitational waves is

$$\begin{equation}\langle {({\Delta}\tilde{{\epsilon}})}^{2}\rangle =\frac{2}{M{\left(\mathfrak{R}\left[{\alpha }_{\text{rel}}^{(1)}\right]+\vert {\beta }_{\text{rel}}\vert \right)}^{2}\left(\mathrm{cosh}\,4{r}_{0}+3\right)}.\end{equation} \tag{ 34 }$$

where r₀ is the initial squeezing of the phonons, to α_rel in the gravitational wave amplitude, with α_rel and β_rel given by equations (29) and (30), respectively. The details of this calculation are provided in appendix A.

In this and all subsequent considerations, we take M = 1. This corresponds to a gravitational wave interacting with a single BEC for N periods. We also assume that the BEC can continuously be regenerated over the course of a year [66, 67], similar to the BEC machine idea suggested by [26], such that observations can be continuously made of a gravitational wave, whose frequency remains approximately constant over the total time of observation. In this case, we find

$$\begin{equation}\sqrt{\langle {({\Delta}\tilde{{\epsilon}})}^{2}\rangle }=\sqrt{\frac{1}{{N}_{\text{tot}}}}\frac{\sqrt{2}}{\left(\mathfrak{R}\left[{\alpha }_{\text{rel}}^{(1)}\right]+\vert {\beta }_{\text{rel}}\vert \right)\sqrt{\mathrm{cosh}\,4{r}_{0}+3}},\end{equation} \tag{ 35 }$$

where N_tot is the total number of regenerations of the condensates. For simplicity, we consider a year of continuous observation, so ${N}_{\text{tot}}\approx \frac{\mathrm{1}\enspace \mathrm{y}\mathrm{r}}{N/f}$.

We note that, in practice, it would be quite difficult to reconstruct the BEC with exactly the same properties. Extra source noise would be introduced into the system, and so such a BEC machine would be beneficial only if the noise resulting from an inexact replication is smaller than the gain in sensitivity resulting from multiple BEC experiments.

We also note that, in equation (35), we are calculating the sensitivity of the scaled gravitational wave amplitude $\tilde{{\epsilon}}$, rather than itself. To obtain ⟨(Δ)²⟩ from $\langle {({\Delta}\tilde{{\epsilon}})}^{2}\rangle$, one can either average or maximize over the components k_x,y of the wavevector. As both methods scale the estimation of the gravitational wave amplitude by $\mathcal{O}(1)$, we will not include this effect in our subsequent calculations and Figures.

As the BEC has a finite size, we note that the minimum gravitational wave frequency that can be observed is $f\approx \frac{2{\bar{c}}_{s}}{L}$ (where we use the relation that resonance between the gravitational wave and phonon occurs at $q=\frac{\omega }{{\Omega}}\approx \frac{1}{2}$ and ${\bar{c}}_{s}$ is the speed of the phonons). The maximum frequency observed is derived from the chemical potential $\mu =m{c}_{s}^{2}$, which sets the upper bound on the frequency of phonons created: f ≪ μ. Again using $q\approx \frac{1}{2}$, we see that the maximum gravitational wave frequency is $f\approx \frac{2m{\bar{c}}_{s}^{2}}{2\pi \hslash \times 10}$, where the 10 comes from the inequality. We show in figure 4 how the frequency bands of our gravitational wave detector are dependent on the speed of sound and length of the condensate. For lower-frequency gravitational wave detection, a larger condensate is required while high-frequency gravitational wave detection requires a faster speed of sound.

Zoom In Zoom Out Reset image size

To illustrate the sensitivity to gravitational waves, we take the number density to be 10²⁰ m⁻³ and a mass of the atoms to be 10⁻²⁵ kg. In figures 5 and 6, we consider the sensitivity of an (undamped) BEC and illustrate the effect of q, length of the condensate, and frequency of incoming gravitational waves on strain sensitivity. We note that a BEC with speed of sound modulation is most sensitive to lower frequency gravitational waves, in contrast to [26, 47]. In section 4, we consider a more realistic situation by modeling Beliaev damping of the phonons within the condensate [68].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This effect of having greater sensitivity away from resonance must be highlighted. When the system is away from $\frac{\omega }{{\Omega}}=\frac{1}{2}\pm \delta$, where δ ≪ 1, we saw that the sensitivity to gravitational waves is actually increased. In these off-resonance cases, the population of phonons grows more slowly, which corresponds to total energy growing more slowly. Hence we can observe gravitational waves for a longer time (a cross-section of this effect is illustrated in figure 5).

We have shown that by using parametric resonance, we can obtain several orders of magnitude in sensitivity (with the level of sensitivity determined by the parameters of the condensate). An additional 3–4 orders of sensitivity can be obtained from observing the gravitational wave over the course of a year. We emphasize that we are implicitly assuming that the gravitational wave's frequency is approximately constant during the total observation time and all the BECs created over the course of a year have the same experimental parameters. If the gravitational wave's frequency changed during the observation time or some of the BECs had slightly different experimental parameters than the other BECs created, this would introduce additional noise sources (beyond what we discuss in section 5) and decrease the sensitivity to gravitational waves.

In section 3, we considered an undamped BEC with no additional sources of noise and demonstrated that, if such a system could be created, then a BEC could potentially be used to detect gravitational waves across several orders of magnitude in frequency, depending on the speed of the phonons and the length of the condensate. In reality, even at zero-temperature (which was implicitly considered), the phonons within the BEC will naturally undergo decoherence through Beliaev damping, where the damping rate is given by [69]

$$\begin{equation}{\gamma }_{B}=\frac{3}{640\pi }\frac{\hslash {\omega }^{5}}{mn{c}_{s}^{5}}\end{equation} \tag{ 36 }$$

for a cubic BEC with no speed of sound modulation, where any correction due to the modulation of the speed of sound would manifest itself as a higher-order term.

Since the system corresponds to a Gaussian state, we can use a covariance matrix to compute the fidelity. As demonstrated in [70], the covariance matrix σ is damped as $\sigma ={\text{e}}^{-{\gamma }_{B}t}{\sigma }_{0}+(1-{\text{e}}^{-{\gamma }_{B}t}){\sigma }_{\infty }$, where σ₀ is the initial covariance matrix and σ_∞ is the t → ∞ covariance matrix. For our purposes, we can neglect the second term as γ_B t ≫ 1 when damping becomes a non-negligible effect. From [63], the quantum Fisher information depends only the square of several combinations of elements of the covariance matrix. Hence, the quantum Fisher information is damped as ${H}_{{\epsilon},\text{damped}}\sim {e}^{-2{\gamma }_{B}t}{H}_{{\epsilon},\text{undamped}}$. Therefore,

$$\begin{equation}\langle {({\Delta}\tilde{{\epsilon}})}^{2}\rangle =\frac{1}{{T}_{\text{tot}}}\frac{2}{M{\left(\mathfrak{R}\left[{\alpha }_{\text{rel}}^{(1)}\right]+\vert {\beta }_{\text{rel}}\vert \right)}^{2}\left(\mathrm{cosh}\,4{r}_{0}+3\right)}{e}^{2{\gamma }_{B}t}\end{equation} \tag{ 37 }$$

where t = N_max T is the maximum time of observation of the condensate, N is the maximum number of oscillations of the speed of sound, and $T=\frac{2\pi }{{\Omega}}$ is the period of the gravitational wave. We note that ${N}_{\mathrm{max}}T< \frac{1}{{\gamma }_{B}}$, where $\frac{1}{{\gamma }_{B}}$ is approximately the decoherence time [70].

In figure 7, we consider how the strain sensitivity is affected when damping is considered for various gravitational wave frequencies and values of q. We see that larger condensates and faster speeds of the phonons give rise to greater sensitivity. However, we note that the sensitivity of a BEC to gravitational waves is optimal at lower frequencies, in contrast to the results of [26, 47]. In figure 7, when fixing q, we see that the strain sensitivity becomes worse further away from the resonance peak of $q=\frac{1}{2}$.

Zoom In Zoom Out Reset image size

4.1. Squeezing

In previous work [47] and in the original proposal of [26], it was noted that initially squeezing the phonon states seemed to be a necessary feature for a BEC to detect gravitational waves. Specifically, to exceed the sensitivity of LIGO, the phonons needed to be squeezed by much more than the capabilities of current experimental designs (7.2 dB), though it was noted that there might be a means [71] of exceeding current limitations.

Though initially squeezing improves the sensitivity to gravitational waves in our current work by a factor of

$$\begin{equation*}\sqrt{3+\mathrm{cosh}\,4{r}_{0}},\end{equation*}$$

we note that this feature is not an essential element of our proposal. For r = 0.83 (potentially achievable with current technology [72]), this improves sensitivity by only a factor of 4.1. Improvements in squeezing may increase the sensitivity to gravitational waves, though we note that squeezed states would contribute to the energy density of the condensate, and therefore cause the condensate to become unstable quicker (section 2.3.2). Furthermore, the squeezing parameter would also decay, as discussed in [70].

At the present time using current techniques (which we discuss in more detail in the next section), it would be difficult to detect gravitational waves using a BEC, though the results of this paper indicate that this could potentially be achieved in the future. When this occurs, then Beliaev damping and the effects of squeezing may need to be considered in more detail. However, in the case of observing a continuous gravitational wave by using parametric resonance, initially squeezing the system and Beliaev damping are not the primary limiting factors as suggested in [47]. Rather it is three-body recombination, which we comment on in the following section.

It is important to note that we have ignored non-zero temperature and quantum depletion effects. As we discuss in the next section, these fluctuations will induce additional sources of noise, which will affect the sensitivity to gravitational waves. In the next section, we briefly consider some of these fluctuations and their effect, though a rigorous analysis is beyond the scope of this work, which we defer to future studies.

5.1. Observable sources

In figure 7, we demonstrated that with a speed of sound of ${\bar{c}}_{s}=5$ cm s⁻¹ and a cubic condensate with side lengths of L = 500 μm, gravitational waves of amplitude $\mathcal{O}(1{0}^{-20})$ could be detected at frequencies around 500 Hz. In figure 5, lower frequency gravitational waves could be detected if their amplitude was $\mathcal{O}(1{0}^{-21})$ and higher. However, as was noted in [25], the amplitude of continuous gravitational waves tends to be weaker than that of transient gravitational waves.

For this reason, it is unlikely that a BEC could detect gravitational waves without going to larger sizes or faster speeds of sound. In figure 8, we illustrate the speeds of sound and condensate lengths necessary to detect gravitational waves with a smaller amplitude, where we also account for Beliaev damping. Currently, it is not possible to create condensates that large, nor is it possible for the condensates to have a long enough lifetime; as we discuss in the next section, larger speeds of sound increase the three-body recombination rate, which in turn decreases the lifetime of the condensate.

Zoom In Zoom Out Reset image size

A BEC is an advantageous gravitational wave detector as it is capable, in principle, of observing gravitational wave sources across several orders of magnitude in frequency, depending on the condensate's length and speed of sound. At gravitational wave frequencies in the tens or hundreds of Hertz, continuous signals from rotating neutron stars or pulsars could be detected [23, 24]. The amplitude of such waves may be $\mathcal{O}(1{0}^{-24})$.

Around black holes, it has been theorized that a boson cloud could form through superradiance, and then emit continuous gravitational waves during annihilation processes [73]. The frequency of the emitted waves will be dependent on the mass of the bosons and distance from Earth, though could potentially be in the kilohertz regime (smaller boson masses will have frequencies in the tens or hundreds of Hertz). The amplitudes of such waves could be as large as $\mathcal{O}(1{0}^{-22})$ [73]. One candidate for boson clouds is axions. In addition to annihilations, axions can undergo energy level transitions, which will also emit gravitational waves. These tend to be in the low-frequency regime ($< 200$ Hz) and have amplitudes $\mathcal{O}(1{0}^{-24})$ and below, though one energy level transition could be $\mathcal{O}(1{0}^{-22})$ (frequency of gravitational wave $\sim 40$ Hz). At higher frequencies, there may be additional cosmological sources, such as cosmic strings or inflationary signatures. However, we note that the signals from such sources may be too weak to be detected by our proposed setup [74].

In figure 8, we also illustrate the amplitude and frequency regimes of such signals. We find that, though it may be necessary to go to larger length scales and speeds of sound to observe most of the continuous gravitational wave sources, it may still be possible to detect gravitational waves using condensates of lengths L ∼ 5 × 10⁻⁴ m and speeds of sound ${\bar{c}}_{s}\sim 5\times 1{0}^{-2}$ m s⁻¹.

5.2. Experimental design

5.3. Experimental challenges

5.3.1. Vibrational noise

In the previous sections, we assumed a cubic trap for simplicity. However, our findings transfer to other trap geometries, as gravitational wave detection is dependent on the existence of particular phonon modes. As long as the desired phonon mode exists, gravitational waves can be detected, regardless of the trap geometry. Different trap geometries will change the quantum Fisher information and the dissipation rate, though the greatest influence will most likely be the trap's size, rather than its specific geometric features.

The dominant source of trap length noise would come from laser intensity noise. Any realistic potential has the shape of a high-order polynomial rather than a perfectly sharp edge. Intentionally modulating the trap laser power in a box trap has been used to excite sound waves in Fermi gases [78]. For this reason, fluctuations in the power of the trapping laser constitute a principle source of noise.

The length of a box trap, in the Thomas–Fermi approximation, is the distance between the points at the two ends where the trapping potential equals the chemical potential. This can be modulated by changing the positions of the light sheets that form the hard walls at either end of the box or by modulating the intensity of the trapping light. To analyse the impact of light intensity noise, we consider the equation $\mu ={U}_{0}{e}^{-2{(z/{w}_{0})}^{2}}$, where μ is the chemical potential, U₀ is the peak potential for a blue-detuned (i.e., light frequency greater than the atomic resonance) light sheet, w₀ is the Gaussian beam waist, and z is the distance from the intensity peak to the edge of the trap defined by the equation. Since U₀ is proportional to the power P in the light sheet, we note that $\mathrm{d}z/\mathrm{d}{U}_{0}={w}_{0}^{2}/(4z{U}_{0})$ implies the change in z, δz, due to a small change in power, δP, is

$$\begin{equation}\delta z=\frac{{w}_{0}^{2}}{4z}\frac{\delta P}{P}.\end{equation} \tag{ 38 }$$

For reasonable laser powers and focal parameters, z ≈ 2w₀. For a cubic box of length L, the minimum waist is roughly found by setting the Rayleigh range $({z}_{\text{R}}=\pi {w}_{0}^{2}/\lambda )$ equal to the length, which gives $\delta z=\sqrt{L\lambda }/(8\sqrt{\pi })\cdot \delta P/P$. Assuming correlated power fluctuations for the sheets at both ends of the trap (as from laser power noise; fluctuations due to beam splitting to make the two sheets should be anti-correlated and so cancel to leading order) we find finally

$$\begin{equation}\frac{\delta L}{L}=\sqrt{\frac{\lambda }{16\pi L}}\frac{\delta P}{P}.\end{equation} \tag{ 39 }$$

Using 532 nm light as in [76] and a length of 1 mm gives a sensitivity δL/L = 0.003(δP/P). The current state of the art for this wavelength could allow for 50 W beams. In an experiment lasting 3 s, assuming the power is controlled to the shot noise limit, we find the length noise is δL/L = 1.5 × 10⁻¹³. This sets a bound on the strain sensitivity for any technique using BEC phonons to detect gravitational waves in a single experiment. To approach the sensitivities needed for actually detecting gravitational waves, different geometries and longer BEC lifetimes would be necessary. We mention some thoughts on bridging these many orders of magnitude below in the future prospects. Note, as mentioned above, that the main conclusions regarding parametric resonance techniques are not tied to the cubic geometry considered in this paper, so the advantages of parametric resonance can be carried forward into more realistic experimental setups in the future.

Mechanical vibrations of the optical setup are the other potential source of trap length noise. To minimize the impact of vibrations, we assume the experiment is performed in a large ultrahigh vacuum chamber with optics built directly into a single piece of Zerodur glass [79] suspended on a vibration control system with comparable specifications as the LIGO mirror suspension system [77]. A phonon mode excited by gravitational waves would have different symmetry than the phonon modes that may be excited by centre of mass oscillations of the trap. However, the frequencies of interest will dominantly transfer as centre of mass oscillations. The strong mismatch between internal vibrational excitation resonances and the pendular modes of the optical system as a whole mean the transfer function from vibrations transmitted through the suspension system to length fluctuations will be much less than one. Assuming the transfer function is 1, we find a strain sensitivity of 10⁻¹⁷. Given that the transfer function can be designed to be well below 1, we find that, even with current technology, vibrational noise can be suppressed below the limits imposed by laser intensity fluctuations. It seems quite feasible to reduce vibrational noise below the level needed to detect gravitational waves with a BEC, thanks to the tremendous advances in this area made by LIGO.

5.3.2. Three-body recombination

Another key challenge, hinted at in the laser noise calculation, is the relatively short lifetime of a BEC due to three-body recombination loss, which has also been discussed in [48]. This is the dominant loss mechanism in most BEC experiments. Three-body recombination occurs when three atoms collide, allowing two to form a molecular bound state and releasing the molecule's binding energy as kinetic energy for the molecule and third atom. This energy will be vastly larger than the trap depth, and so both atom and molecule are ejected from the trap and lost.

We use 3 s as a standard for lifetime of a BEC due to three-body loss, because for alkali atoms at roughly the densities called for (set by the desired speed of sound), lifetimes would be expected to range from about one to ten seconds. However, three-body loss does not lead to a simple exponential decay as implied by assigning a lifetime. Rather, since it is a three body process, we have $\dot{n}=-{K}_{3}{n}^{3}$, which leads to $n(t)={n}_{0}/\sqrt{{n}_{0}^{2}{K}_{3}t+1}$, where n₀ = n(0), assuming no other loss processes are substantial, where ${K}_{3}$ is the three-body loss coefficient.

The rapid loss early in experiments would cause the speed of sound to change on the time scales relevant to gravitational wave modulations. The changing speed of sound would also change the phonon dispersion relation, destroying the coherence of the gravitational wave coupling. This could drop the effective length of an experiment to well below 1 s, since all BEC-phonon schemes for gravitational wave detection rely on coherent coupling of the gravitational wave to the phonon modes over long time periods. On time scales of 1–10 s, we can possibly solve this problem through a tailored ramp of the interaction strength. Since ${c}_{s}=\sqrt{gn/m}$, we can compensate a known form for density loss over time by increasing the interaction strength, g, in a complementary fashion to keep c_s fixed. This ramp of the interaction strength can be accomplished by tuning the ambient magnetic field through the same physical mechanism by which we modulate the interaction strength to produce parametric resonance.

While this sweep of g can fix the initial problem of three-body loss limiting coherence times to well below 1 s, it cannot be used for arbitrarily long times. Since precision scales as N^−1/2, the loss of atoms reduces precision. Further, increasing g increases K₃. Roughly speaking, we can expect K₃ ∝ g⁴, so if we ramp g to hold gn fixed, the quantity K₃ n³ will increase like 1/n, accelerating the atom loss. Accounting for all of this, we see 3 s is a good general purpose estimate for the limits of coherence time and lifetime of a BEC for detecting gravitational waves. In future prospects we mention one idea for extending this lifetime that goes beyond the scope of this paper.

In figure 9, we depict the sensitivity curve of a BEC to a 1000 Hz gravitational wave, where the observation time is taken to be the minimum of the non-linearity time and the 3 s lifetime of the BEC. The minimum of the sensitivity curve corresponds to the location where the non-linearity time is exactly equal to the lifetime of the BEC. Moving outside the resonance bands causes the non-linearity time to exceed the lifetime of the condensate, which implies that there is insufficient time to use parametric resonance to amplify the sensitivity. In this case, such a BEC would be insensitive to gravitational waves.

Zoom In Zoom Out Reset image size

5.3.3. Techniques for improvement and other sources of noise

The trap length noise can potentially be improved by squeezing the intensity noise of the laser, since the phase of the trapping light is not critical. State of the art squeezing achieves 15 dB reduction in noise [80, 81]. Supposing a substantial improvement to 30 dB is possible in the future, the achievable strain sensitivity of $\sim 1{0}^{-16}$ is still far behind the sensitivity of LIGO. Using an asymmetric box potential and only looking at phonon modes along the long direction would allow for larger L and a smaller beam waist for the walls in the important direction. A 5 mm long trap with transverse width of 100 μm could improve the sensitivity by nearly an order of magnitude. We conclude that the limit imposed by plausible near-term technologies is a strain sensitivity of 10⁻¹⁷ per shot.

Our calculations in the preceding sections have assumed a zero-temperature BEC. This is a simplifying approximation, though we do not expect temperature-dependence to yield significantly different results. As discussed in [82], a temperature of 150 nK yielded comparable sensitivities as a 0 nK BEC. We also note that BEC temperatures can be on the scale of nK or pK, which would therefore minimize temperature-dependent noise effects [83]. Perhaps more serious is the fact that by parametrically driving the condensate, there will be additional heating effects, which will, over time, change the speed of sound. However the key to our proposal is not the speed of sound itself but the frequency of the speed of sound fluctuations. As long as the speed of sound fluctuations are on resonance at the gravitational wave frequency, a BEC could still be used to detect gravitational waves, though the sensitivity to the waves would change. As illustrated in equation (31), ${\bar{c}}_{s}$ helps set the limit to the maximum number of oscillations N, and thus the maximum parametric resonance time of the BEC. Therefore, this heating effect could provide stronger limits to the maximum time that the BEC can undergo parametric resonance. We leave detailed calculations of this source of noise for future investigation.

We have discussed the noise in the trap, though additional noise will be introduced in the measurement portion of such an experiment. Even if the technical noise of the trap is mitigated, there will remain a few caveats. In section 4, we noted that faster speeds of sound gave rise to greater sensitivity. This was because the damping rate was ${\gamma }_{B}\propto {c}_{s}^{-5}$. With faster speeds of sounds, a BEC will have a larger minimum observable gravitational wave frequency (∝c_s ), though a much larger maximum observable gravitational wave frequency $(\propto {c}_{s}^{2})$. However, at larger speeds of sound, we note that the maximum time of observation will no longer be given by the non-linearity condition of section 2.3.2, but by three-body recombination, as we just discussed. Larger speeds of sound correspond to greater number densities, implying a greater decrease in the number atoms in the condensate [84, 85]. We also note quantum depletion, which describes the loss of atoms from the condensate because of interactions between the atoms and excitations, will contribute to additional noise in the condensate, where the number of atoms in an excited state (therefore, leaving the condensate) is given by [86]

$$\begin{equation}\frac{{n}_{\text{ex}}}{N}=\frac{8}{3\sqrt{\pi }}{(n{a}^{3})}^{1/2}\end{equation} \tag{ 40 }$$

where a is the scattering length and is valid for a zero-temperature uniform Bose gas.

We have considered a BEC with oscillating speed of sound to investigate its sensitivity toward a continuous gravitational wave of the form h₊ =  sin Ωt. We showed that the oscillation gives rise to parametric resonance in the system, which can be exploited to improve the sensitivity to gravitational waves at lower frequencies compared to higher frequencies.

Despite the attractive aspects of this idealized setup, this proposal has formidable experimental challenges. The technical noise present using current technology would overwhelm signals from likely gravitational wave sources. The finite lifetime of real BECs also imposes a difficult constraint, which we implemented as a hard limit of 3 s on the length of an experiment. However, if the technical or lifetime issues can be overcome, our results indicate that parametric resonance can make a BEC a more promising gravitational wave detector than earlier considerations have indicated [44, 47] For simplicity, we have considered a cubic condensate. In more realistic trap geometries, the same analysis can be applied once it is determined which phonon mode is needed for gravitational wave detection.

We offer three prospects for improving the sensitivity beyond the limitations discussed on the paper. Three-body loss coefficients are relatively similar among laser coolable atoms. The one atom that has produced a BEC that could potentially have a dramatically longer lifetime is hydrogen. Since the triplet potential of H₂ supports no bound states, there is no three-body loss in spin-polarized H. Hydrogen BEC experiments in magnetic traps suffer from dipolar relaxation [87]. In principle, H polarized in the lowest-energy Zeeman sublevel could produce BECs with very long lifetimes. This state could be trapped in an optical potential, though using repulsive light sheets as described in the paper would be infeasible, given the unavailability of strong laser sources at wavelengths shorter than 121 nm. Using an attractive potential produced by a CO₂ laser could lead to lifetimes well over a minute.

State of the art CO₂ lasers can produce more than 10 kW of power. Using 25 kW of 10.6 μm light produced by a CO₂ laser would give a shot-noise limit 100 times lower than the 50 W of 532 nm light discussed in the paper. Analyzing the impact of laser intensity noise in an attractive, harmonic trap is beyond the scope of this paper. However, even the repulsive optical potentials described in this paper could potentially have far lower δP/P by using power build-up cavities to produce the trapping potentials.

Finally, at the price and size scale of present gravitational detectors, one could imagine a set of machines producing many BECs for simultaneous measurements to further improve sensitivity. Though using a BEC to detect gravitational waves will be extremely difficult, it should be noted that LIGO was first conceptualized in the 1980s and only detected gravitational waves in 2015. BECs have a head start in knowing that gravitational waves can be directly detected and some of LIGO's technical innovations could potentially be adapted for use by a BEC GW detector. Though such an experiment may take years or decades to achieve, doing so will provide a new method to detecting gravitational waves across potentially several orders of magnitude in frequency.

We thank Aditya Dhumuntarao, Soham Mukherjee, and Erickson Tjoa for their useful comments and discussions. We also thank the participants of the mini-symposium on Quantum Sensors and Their Applications in Fundamental Physics and Astrophysics Experiments (December 13, 2019) at Washington University in St. Louis for their helpful questions and suggestions. We thank Vladimir Dergachev for his valuable discussion regarding LIGO's calculations. We also thank Denis Martynov for providing LIGO's noise spectrum. MR was funded by a National Science and Engineering Research Council of Canada (NSERC) graduate scholarship. AOJ was funded by the University of Waterloo and the Institute for Quantum Computing. This research was supported in part by NSERC and the Perimeter Institute for Theoretical Physics, and by AOARD Grant FA2386-19-1-4077. Research at Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Economic Development & Innovation.

No new data were created or analysed in this study.

As mentioned in section 3, an unknown parameter in a quantum system can be estimated with the quantum Cramer–Rao bound [29],

$$\begin{equation}\langle {({\Delta}\tilde{{\epsilon}})}^{2}\rangle \geqslant \frac{1}{M{H}_{\tilde{{\epsilon}}}},\end{equation} \tag{ 41 }$$

with $\langle {({\Delta}\tilde{{\epsilon}})}^{2}\rangle$ the expectation value of the uncertainty in $\tilde{{\epsilon}}$, M the number of measurements of the system, and ${H}_{\tilde{{\epsilon}}}$ the Fisher information [29]

$$\begin{equation}{H}_{\tilde{{\epsilon}}}=\frac{8\left(1-\sqrt{F({\rho }_{\tilde{{\epsilon}}},{\rho }_{\tilde{{\epsilon}}+d\tilde{{\epsilon}}})}\right)}{d{\tilde{{\epsilon}}}^{2}},\end{equation} \tag{ 42 }$$

where $F({\rho }_{\tilde{{\epsilon}}},{\rho }_{\tilde{{\epsilon}}+d\tilde{{\epsilon}}})={\left[\mathrm{T}\mathrm{r}\,\sqrt{{\rho }_{\tilde{{\epsilon}}}\sqrt{{\rho }_{\tilde{{\epsilon}}+d\tilde{{\epsilon}}}}{\rho }_{\tilde{{\epsilon}}}}\right]}^{2}$ is the fidelity [64, 65]. Note that we can recast equations (29) and (30) as ${\alpha }_{\text{rel}}={\text{e}}^{-\text{i}{\theta }_{\alpha }}\mathrm{cosh}\,{r}_{\text{rel}}$ and ${\beta }_{\text{rel}}={\text{e}}^{-\text{i}{\theta }_{\beta }}\mathrm{sinh}\,{r}_{\text{rel}}$, where r_rel is the squeezing that results from the gravitational wave and ${\text{e}}^{-\text{i}{\theta }_{\alpha }},{\text{e}}^{-\text{i}{\theta }_{\beta }}$ are phase factors. Therefore, we find ${r}_{\text{rel}}=\mathrm{log}\left[\vert {\alpha }_{\text{rel}}\vert +\vert {\beta }_{\text{rel}}\vert \right]$.

Let us suppose that before a gravitational wave interacts with the condensate, our phonons are in a squeezed state, with squeezing parameter r₀ and quadrature angle ϕ₀, where the state of the phonons is $\vert {\zeta }_{0}\rangle ={S}_{0}\vert 0\rangle =\mathrm{exp}\left[\frac{1}{2}({\zeta }_{0}^{\ast }{\hat{a}}^{2}-{\zeta }_{0}{\hat{a}}^{{\dagger}2})\right]\vert 0\rangle$, S₀ is the squeezing operator, and ${\zeta }_{0}={r}_{0}{\text{e}}^{\text{i}{\phi }_{0}}$

Let ${{\epsilon}}_{1}=\tilde{{\epsilon}}$ and ${{\epsilon}}_{2}=\tilde{{\epsilon}}+d\tilde{{\epsilon}}$. For a pure state, the quantum Fisher information is ${H}_{\tilde{{\epsilon}}}=\frac{8(1-\vert \langle {{\epsilon}}_{1}\vert {{\epsilon}}_{2}\rangle \vert )}{{({{\epsilon}}_{1}-{{\epsilon}}_{2})}^{2}}$. The state after interaction with a gravitational wave is $\vert {{\epsilon}}_{i}\rangle ={S}_{{{\epsilon}}_{i}}{S}_{0}\vert 0\rangle$, where ${S}_{{{\epsilon}}_{i}}$ encodes the gravitational wave's influence. Then, ${S}_{\text{rel}}{:=}{S}_{{{\epsilon}}_{1}}^{{\dagger}}{S}_{{{\epsilon}}_{2}}$, so that

$$\begin{equation}\langle {{\epsilon}}_{1}\vert {{\epsilon}}_{2}\rangle =\langle 0\vert {S}_{0}^{{\dagger}}{S}_{\text{rel}}{S}_{0}\vert 0\rangle \end{equation} \tag{ 43 }$$

$$\begin{equation}=\langle 0\vert {S}_{0}^{{\dagger}}\,\mathrm{exp}\left[\frac{1}{2}({\zeta }_{\text{rel}}^{\ast }{\hat{a}}^{2}-{\zeta }_{\text{rel}}{\hat{a}}^{{\dagger}2})\right]{S}_{0}\vert 0\rangle \end{equation} \tag{ 44 }$$

where ${\zeta }_{\text{rel}}={r}_{\text{rel}}{\text{e}}^{\text{i}{\phi }_{\text{rel}}}$. Noting that ${r}_{\text{rel}}\propto \tilde{{\epsilon}}$, we have

$$\begin{align}\hfill \langle {{\epsilon}}_{1}\vert {{\epsilon}}_{2}\rangle & \sim 1+\frac{1}{2}\langle 0\vert {S}_{0}^{{\dagger}}\left({\zeta }_{\text{rel}}^{\ast }{\hat{a}}^{2}-{\zeta }_{\text{rel}}{\hat{a}}^{{\dagger}2}\right){S}_{0}\vert 0\rangle +\frac{1}{8}\langle 0\vert {S}_{0}^{{\dagger}}{\left({\zeta }_{\text{rel}}^{\ast }{\hat{a}}^{2}-{\zeta }_{\text{rel}}{\hat{a}}^{{\dagger}2}\right)}^{2}{S}_{0}\vert 0\rangle .\hfill \end{align} \tag{ 45 }$$

Using ${S}_{0}^{{\dagger}}a{S}_{0}=a\,\mathrm{cosh}\,{r}_{0}-{a}^{{\dagger}}{\text{e}}^{\text{i}{\phi }_{0}}\,\mathrm{sinh}\,{r}_{0}$ and ${S}_{0}^{{\dagger}}{a}^{{\dagger}}{S}_{0}={a}^{{\dagger}}\,\mathrm{cosh}\,{r}_{0}-a{\text{e}}^{-\text{i}{\phi }_{0}}\,\mathrm{sinh}\,{r}_{0}$, we find

$$\begin{equation}\langle 0\vert {S}_{0}^{{\dagger}}{a}^{{\dagger}2}{S}_{0}\vert 0\rangle =-{\text{e}}^{-\text{i}{\phi }_{0}}\,\mathrm{cosh}\,{r}_{0}\,\mathrm{sinh}\,{r}_{0}\end{equation} \tag{ 46 }$$

$$\begin{equation}\langle 0\vert {S}_{0}^{{\dagger}}{a}^{2}{S}_{0}\vert 0\rangle =-{\text{e}}^{\text{i}{\phi }_{0}}\,\mathrm{cosh}\,{r}_{0}\,\mathrm{sinh}\,{r}_{0}\end{equation} \tag{ 47 }$$

$$\begin{equation}\langle 0\vert {S}_{0}^{{\dagger}}{a}^{{\dagger}4}{S}_{0}\vert 0\rangle =3{e}^{-2i{\phi }_{0}}\,{\mathrm{cosh}}^{2}\,{r}_{0}\,{\mathrm{sinh}}^{2}\,{r}_{0}\end{equation} \tag{ 48 }$$

$$\begin{equation}\langle 0\vert {S}_{0}^{{\dagger}}{a}^{4}{S}_{0}\vert 0\rangle =3{e}^{2i{\phi }_{0}}\,{\mathrm{cosh}}^{2}\,{r}_{0}\,{\mathrm{sinh}}^{2}\,{r}_{0}\end{equation} \tag{ 49 }$$

$$\begin{equation}\langle 0\vert {S}_{0}^{{\dagger}}{a}^{2}{a}^{{\dagger}2}{S}_{0}\vert 0\rangle =2\,{\mathrm{cosh}}^{4}\,{r}_{0}+{\mathrm{cosh}}^{2}\,{r}_{0}\,{\mathrm{sinh}}^{2}\,{r}_{0}\end{equation} \tag{ 50 }$$

$$\begin{equation}\langle 0\vert {S}_{0}^{{\dagger}}{a}^{{\dagger}2}{a}^{2}{S}_{0}\vert 0\rangle =2\,{\mathrm{sinh}}^{4}\,{r}_{0}+{\mathrm{cosh}}^{2}\,{r}_{0}\,{\mathrm{sinh}}^{2}\,{r}_{0}.\end{equation} \tag{ 51 }$$

Therefore,

$$\begin{equation}\vert \langle {{\epsilon}}_{1}\vert {{\epsilon}}_{2}\rangle \vert \sim 1-\frac{3{r}_{\text{rel}}^{2}}{16}-\frac{1}{16}{r}_{\text{rel}}^{2}\,\mathrm{cosh}(4\text{r0})+\frac{1}{8}{r}_{\text{rel}}^{2}\,\mathrm{cos}(2{\phi }_{0}-2{\phi }_{\text{rel}}){\mathrm{sinh}}^{2}(2{r}_{0}).\end{equation} \tag{ 52 }$$

Therefore, after averaging over angles, we find

$$\begin{equation}{({{\epsilon}}_{1}-{{\epsilon}}_{2})}^{2}{H}_{{\epsilon}}=\frac{1}{2}{r}_{\text{rel}}^{2}\left(3+\mathrm{cosh}({r}_{0})\right).\end{equation} \tag{ 53 }$$

Recalling that ${r}_{\text{rel}}=\mathrm{log}\left[\vert {\alpha }_{\text{rel}}\vert +\vert {\beta }_{\text{rel}}\vert \right]$ and from equations (29) and (30),

$$\begin{equation}\ {\alpha }_{\text{rel}}^{\ast }=1+\tilde{{\epsilon}}\left({\alpha }_{\tilde{{\epsilon}}}^{\ast }{\alpha }_{2a}-{\beta }_{\tilde{{\epsilon}}}^{\ast }{\beta }_{2a}\right),\end{equation} \tag{ 54 }$$

$$\begin{equation}{\beta }_{\text{rel}}^{\ast }=-\tilde{{\epsilon}}\left({\alpha }_{\tilde{{\epsilon}}}^{\ast }{\beta }_{2a}^{\ast }-{\beta }_{\tilde{{\epsilon}}}^{\ast }{\alpha }_{2a}^{\ast }\right),\end{equation} \tag{ 55 }$$

we find ${r}_{\text{rel}}\sim {\left(\mathfrak{R}\left[{\alpha }_{\text{rel}}^{(1)}\right]+\vert {\beta }_{\text{rel}}\vert \right)}^{2}{{\epsilon}}^{2}$, where ${\alpha }_{\text{rel}}=1+{\epsilon}{\alpha }_{\text{rel}}^{(1)}$ Therefore, with equations (41) and (53), the sensitivity to gravitational waves is given by

$$\begin{equation}\langle {({\Delta}\tilde{{\epsilon}})}^{2}\rangle =\frac{2}{M{\left(\mathfrak{R}\left[{\alpha }_{\text{rel}}^{(1)}\right]+\vert {\beta }_{\text{rel}}\vert \right)}^{2}\left(3+\mathrm{cosh}(4{r}_{0})\right)}.\end{equation} \tag{ 56 }$$