We introduce a quantum interferometric scheme that uses states that are sharp in frequency and delocalized in position. The states are frequency modes of a quantum field that is trapped at all times in a finite volume potential, such as a small box potential. This allows for significant miniaturization of interferometric devices. Since the modes are in contact at all times, it is possible to estimate physical parameters of global multimode channels. As an example, we introduce a three-mode scheme and calculate precision bounds in the estimation of parameters of two-mode Gaussian channels. This scheme can be implemented in several systems, including superconducting circuits, cavity-QED, and cold atoms. We consider a concrete implementation using the ground state and two phononic modes of a trapped Bose–Einstein condensate. We apply this to show that frequency interferometry can improve the sensitivity of phononic gravitational waves detectors by several orders of magnitude, even in the case that squeezing is much smaller than assumed previously, and that the system suffers from short phononic lifetimes. Other applications range from magnetometry, gravimetry, and gradiometry to dark matter/energy searches.

## I. INTRODUCTION

Interferometers have become a powerful tool for precision measurements, often achieving sensitivities that are not possible using any other known technique. In the most common implementation, waves travel along two different spatial paths and are recombined creating an interference pattern. We will call this setup a *spatial interferometer* since the system follows two different trajectories in space. The phase difference along the paths is a *U*(1) channel that encodes physical parameters, such as frequency and field strengths, that can be measured with very high precision. Interferometers can operate in both the classical and quantum regimes. The most remarkable application of a spatial interferometer is perhaps LIGO, which uses a very similar interferometer to that designed by Michelson to observe gravitational waves.^{1} Early measurements used classical light, but recently, the sensitivity has been enhanced by injecting quantum states of light through the output Faraday isolator.^{2}

Interferometry in the quantum regime uses not only photons but also the wave character of massive particles, such as electrons, neutrons, atoms, and molecules. Here, each particle generically follows a superposition of two spatial trajectories and interferes with itself when the paths recombine. The phase difference acquired by a spatial quantum superposition has been used, for example, to measure accelerations;^{3} fundamental constants,^{4} including the fine-structure constant;^{5} and set constraints on dark energy models;^{6} and has also been proposed as a method to detect gravitational waves at low frequencies.^{7} For a review of cold-atom interferometry, see, for example, Refs. 8–11.

In spatial interferometry, cutting-edge sensitivity is usually limited by the available time of flight. This generically requires large spatial separations or long interferometer arms, which can be several meters or kilometers long. Typically, spatial interferometers cannot be reduced in size without losing precision. However, the time of flight in atom interferometry can be increased by using Bragg diffraction and Bloch oscillations to slow down the particles^{12} as well as introducing multiple loops^{13} and guided schemes for rotation sensing.^{14,15} In such schemes, interactions are undesirable because they reduce the coherence time of the interferometer.^{16} Here, we propose using phonons, which appear only in the presence of atom–atom interactions, as a way to miniaturize devices while keeping high precision. The idea is to use interferometry in the frequency domain, or, equivalently, temporal domain, and we refer to this as *quantum frequency interferometry*. In this case, waves do not follow different spatial paths but are quantum modes of vibration, such as those produced by interacting particles trapped in a localized potential. Since the modes are non-local and they can interact at all times, it is also possible not only to estimate quantities encoded in phase channels but also to estimate the parameters of global unitaries, including entangling operations, in the Hilbert space of the modes.

This new type of quantum interferometry is inspired by recent non-interferometric studies that consider estimating physical parameters using phonons of Bose–Einstein condensates (BECs). For example, frequency modes of phonons were considered in theoretical studies for measuring accelerations^{17} and to detect high-frequency gravitational waves.^{18,19} More recently, phonons have been used in theoretical studies that consider miniaturizing gravimeters^{20} and measuring the gravitational field gradient within the millimeter scale.^{21} In this paper, we extend these ideas to frequency interferometry (FI), the application of which can improve by several orders of magnitude the sensitivities reported in these studies, which we illustrate using the phononic gravitational wave detector.^{18}

## II. SPATIAL INTERFEROMETRY AND THE PUMPED-UP *SU*(1, 1) SCHEME

In general, a quantum interferometer can be broken up into four stages: (i) the active/passive beam splitter, which usually entangles/populates the input state, (ii) the channel(s) that imprint(s) the parameter to be estimated, and (iii) the transformation that recombines the modes, facilitating interference. Finally, in (iv), the modes are measured. Conventionally, the second stage is considered to consist of unitaries acting independently in each spatial mode $U\u0302(\varphi )=U1(\varphi )\u2297\u2009U2(\varphi )\u2297\cdots $ (see also Fig. 1).

Spatial interferometers, including the well-known Mach–Zehnder interferometer, are commonly SU(2) interferometers where the states belong to a two-dimensional Hilbert space, which basis is spanned by the states corresponding to the two possible paths $|1\u3009$ and $|2\u3009$ followed by the particles, and all the elements of the interferometer can be described in terms of quantum operators that obey an SU(2) algebra.^{22} For example, beam splitters exchange particles from one path to another and can be mathematically represented by $\sigma +=|1\u3009\u30082|$ and $\sigma \u2212=|2\u3009\u30081|$. In 1986, Yurke *et al.* introduced a new type of interferometer where the passive beam splitters of SU(2), which perform mode-mixing operations, are replaced with active beam splitters, which perform quantum squeezing operations, such that the algebra that defines its operations is SU(1,1).^{22} This new type of interferometry, called SU(1,1) interferometry, has the advantage over SU(2) that squeezing and entanglement, which can improve the sensitivity of an interferometer, are generated within the interferometer itself rather than this having to be sourced separately. Then, since the overall scheme will have fewer operational elements, it should, in principle, be more robust to noise. In this case, squeezing between spatial modes is introduced at points $p1$ and $p2$ in Fig. 1.

Sending in a coherent pump to an SU(1,1) interferometer results in a two-mode squeezed state in the two spatial modes and, therefore, in principle, a $1/N$ scaling in sensitivity compared to $1/N$ for a classical side mode state. However, generating a large number of particles in the spatial modes is extremely challenging, and so the sensitivity is easily beaten by interferometers operating at the standard quantum limit with large input states. In order to overcome this issue, a variant of the SU(1,1) interferometer, called pumped-up SU(1,1), has recently been proposed where the pump beam first goes through a parametric amplifier and is then mixed with the side modes, such that all particles take part in the estimation procedure.^{23} This essentially allows for a $1/NN0$ scaling in sensitivity, where *N*_{0} is the number of particles in the pump beam and, in general, $N0\u226bN$. This interferometer contains both passive and active beam splitters, such that all the beam splitters within it are generated by general two and three-mode unitary Gaussian operations. In some sense, this interferometer is a generalization of SU(2) and SU(1,1) interferometry, with the parametric amplifiers of the SU(1,1) scheme seeding the quantum input state of an SU(2) (or more properly SU(3) since there are three modes) interferometer. Depending on the chosen three-mode-mixing angle, it is also able to act like an SU(1,1) or SU(2) interferometer. After the three-mode-mixing operation, the modes are physically separated and undergo unitary transformations that encode phases $\varphi 1$ and $\varphi 2$ on the respective states. For example, in an SU(2) interferometer, the relative phase $\varphi 1\u2212\varphi 2$ is estimated, whereas an SU(1,1) interferometer is sensitive to the total unitary transformation $U\u0302(\varphi )=exp\u2009(\u2212i\varphi N\u0302/2)$, where $\varphi :=\varphi 1+\varphi 2$ and $N\u0302:=a\u03021\u2020a\u03021+a\u03022\u2020a\u03022$, with $a\u03021$ and $a\u03022$ the annihilation operators for the two side modes^{22} (see Fig. 1). Note that the unitary can be written as two transformations acting independently on each mode $U\u0302(\varphi )=e\u2212i\varphi 2a\u03021\u2020a\u03021e\u2212i\varphi 2a\u03022\u2020a\u03022=U1\u2297U2$.

## III. QUANTUM FREQUENCY INTERFEROMETRY

In matter-wave spatial interferometers, interactions are usually suppressed, so that all particles are independent of each other. In that sense, they are single particle setups. The initial state is a product state of the individual particle states, and the sensitivity scales with the standard quantum limit $1/N$. However, quantum metrology studies show that entanglement, which requires interactions, can be used to increase precision, reaching, in the optimal case, the Heisenberg limit $1/N$.^{24} Interactions are key to frequency interferometry: the atomic interactions give rise to collective excitations, which are described by massless bosons (also known as phonons) that are used as the quantum information carriers. We consider initially *N* frequency phonon modes of a localized potential (Fig. 2, top). First, the modes are prepared, often in the vacuum, and then a passive or active beam splitter $U\u0302(\theta )$, with $\theta =(\theta 1,\theta 2,\u2026,\theta N)$, is applied to create a quantum state of a subset of modes, which usually involves populating and entangling these modes. Subsequently, the modes (or an even smaller subset of modes) go through a channel $U\u0302(\u03f5)$ with $\u03f5=(\epsilon 1,\epsilon 2,\u2026,\epsilon m)$, which imprints the parameters to be estimated. Finally, before measurement, the modes are “recombined” by $U\u0302(\u2212\theta )$. The frequency modes of the phonons are delocalized in the potential but in physical contact with each other, allowing for global estimation channels where a global channel is a unitary that cannot be written as $U(\epsilon )\u2260U1\u2297U2\u2297\u2026$. This, in particular, facilitates the generalization of pumped-up *SU*(1, 1) interferometry to global channels that include unitaries that entangle the modes. Sensitivities are increased by preparing entangled states, which can be achieved using the particle interactions. The precision depends on the lifetime of the mode excitations, which can be extended by tuning the interactions without making the interferometer larger.

### A. Frequency interferometry with Gaussian channels

Here, we have chosen to present a three-mode example of frequency interferometry, illustrated in Fig. 3. The example uses an analogue in frequency space of the squeezing and tritter operations used in the pumped-up SU(1,1) scheme. However, the pumped-up SU(1,1) was introduced for separable channels of the form $U\u0302(\varphi )=U1\u2297U2$. In this example, the channels we consider are entangling two-mode Gaussian channels given by two-mode squeezing and mode-mixing channels,

or

respectively, where $\xi :=sei\varphi B$ and $\zeta :=mei\varphi A$, with $s\u22650,\u2009m\u22650$, and $\varphi A,\varphi B\u2208\mathbb{R}$. These channels include both SU(1,1) and SU(2)-like interferometry. Together with the phase-shift channel, considered in Ref. 23 and single-mode squeezing, they form the complete set of unitary two-mode Gaussian channels, also known as a Bogoliubov transformations.^{25} An analysis of the optimum Gaussian input states for Gaussian channels using the Quantum Fisher Information (QFI) is given in Ref. 26.

We will consider the parameter of interest *ε* to be encoded in the squeezing parameter *s* and mode-mixing parameter *m*. Specifically, we take $s=:14\epsilon B$ and $m=:14\epsilon A$, where *A* and *B* depend on the physical quantities of the specific implementation of the scheme. Since the input is a Gaussian state and all operations within the interferometer are Gaussian, we find it most straightforward to consider the operation of the new interferometer within the covariance matrix formalism (CMF) (see, e.g., Ref. 25), where just two finite matrices are needed to describe the quantum states of the modes. This formalism allows for simple calculations of the QFI for the device, where we surprisingly find that the sensitivity can be slightly enhanced when the parameter of interest is encoded in a squeezing or mode-mixing channel.

Although the QFI can provide the optimum precision given the above first two stages of an interferometer, it does not by itself identify the measurement scheme that achieves this. However, we find that the sensitivity of a simple intensity measurement scheme can saturate the quantum Cramér–Rao bound (QCRB) at large input particle number. Within the CMF, we also find simple expressions for the sensitivity of general interferometers with such intensity measurement schemes, which can be straightforwardly generalized to other schemes such as homodyne.

The CMF is a phase-space representation of a quantum state where a Gaussian state is fully defined by its displacement vector *d* and covariance matrix $\sigma $. In the real *q*–*p* representation, these are defined as the following for a system consisting of *n* bosonic modes:^{32}

where $x\u0302:=(x\u03021,x\u03022,\u2026,x\u03022n\u22121,x\u03022n)T$ and $x\u0302i$ are quadratures defined by

where $i\u2208\u2009Z+$ and $a\u0302i$ and $a\u0302i\u2020$ are the annihilation and creation operators, respectively. Unitary transformations ** U** acting on density matrices now lead to symplectic matrices

**acting on the displacement and covariance matrices through $d\u2032=Sd$ and $\sigma \u2032=S\sigma ST$.**

*S*^{25}

The initial state of the pump mode is assumed to be a coherent state. As in standard SU(1,1) interferometry, we act on this state with a two-mode squeezing operation $U\u0302sq(r)=exp\u2009{\chi (a\u03021\u2020a\u03022\u2020\u2212a\u03021a\u03022)}$ to parametrically populate the side modes, where $\chi :=r\u2009exp\u2009{i\u03d1sq}$. The state of the full system is then given by $Ssd0$ and $Ss\sigma 0SsT$, where $Ss$ is the symplectic matrix of this squeezing unitary and $d0$ and $\sigma 0$ are the displacement and covariance matrices of the initial coherent state, respectively.^{33} Here, we have assumed that the pump is fairly undepleted by the squeezing operation and remains in a coherent state $|\alpha \u3009$, but we take $\alpha \u2192\alpha 0$ after acting with $Ss$, where $|\alpha 02|:=|\alpha |2\u22122\u2009sinh2\u2009r$ and $N0:=|\alpha 0|2,\u2009N:=2\u2009sinh2\u2009r$, so that particle number is conserved.^{23}

Next, we apply a tritter to the three modes, whose symplectic matrix we denote as $Str$. We then act on the side modes with the squeezing or mode-mixing operations given by (1) and (2). Subsequently, the beams are brought “back together” with another tritter and then an outcoupling process, which are both the reverse of the operations that were performed prior to the Gaussian unitary channel. The state of the full interferometer is then defined by $d=Sd0$ and $\sigma =S\sigma 0ST$, where $S:=S\u2212S\epsilon S+$ with $S\u2212:=Ss(\u2212r)St(\u2212\theta ),S+:=St(\theta )Ss(r)$, and $S\epsilon $ being either the squeezing or mode-mixing channel for the side modes.

### B. Quantum Fisher information

Since it is independent of the particular measurement scheme used, when calculating the QFI, we only need to consider the operations up to and including the Gaussian unitary channels, i.e., the state of the relevant system is defined by $d=S\epsilon S+d0$ and $\sigma =S\epsilon S+d0(S\epsilon S+)T$. For Gaussian states, the QFI and $H\epsilon $ can be obtained through simple expressions (see, e.g., Refs. 34 and 35). If the parameter of interest is encoded through a squeezing channel, the QFI is

where

In Fig. 4, we plot the dependence on the tritter angle *θ* for the QFI when there is a squeezing channel (7), standard phase-shift channel (which can be found in Ref. 23 and Appendix C), and mode-mixing channel (8). Interestingly, we find a slight improvement in the QFI for the squeezing and mode-mixing channels compared to the conventional phase-shift channel. In the former case, this is because even a vacuum input to the channel can be used to estimate the parameter of interest.

Given the optimum phase relationships $\u03d1sq=\varphi B+\pi /2$ and $\nu B=0$, the QFI for the squeezing channel (7) has turning points at *θ* = 0, $\theta \u2212\pi /2$, and $\theta =\theta t$, where *θ _{t}* is defined in Ref. 23 and Appendix D and approximates $\pi /4+csc\u22121(N+N(N+2))$ when $N\xaf$ is large. Here, we concentrate on the most interesting and relevant regimes, which are when

*θ*= 0 and $N\xaf\u226b1$. In the former case, we recover standard SU(1,1) interferometry, and so (7) becomes the QFI for an SU(1,1) interferometer with a squeezing channel, which could still be considered an SU(1,1) interferometer since the unitary representation of a squeezing channel is part of the SU(1,1) group, and the number-sum operation can still be used in the measurement scheme. The QFI in this case scales as

*N*

^{2}, which is the scaling obtained in conventional SU(1,1) interferometry. On the other hand, when $N\xaf\u226b1$, the QFI (7) can be approximated by $B2\u2009sin2(2\theta )NN\xaf/8$, where we have also assumed that $N\u226b2$ and taken the optimum phase relation $\nu B=0$. This $NN\xaf$ scaling was also found in Ref. 23 for the phase-shift channel case. In practice, this scaling can beat the conventional

*N*

^{2}scaling by orders of magnitude since $N\xaf\u226bN$ in current experiments.

If instead of the squeezing channel, the parameter of interest is encoded through a mode-mixing channel, the QFI is

where

As for the squeezing channel, when *θ* = 0, we obtain the QFI for a conventional SU(1,1) interferometer with a mode-mixing channel. In contrast to the squeezing case, however, this type of interferometer derived when *θ* = 0 would not be considered an SU(1,1) interferometer by the original definition^{22} since the unitary representation of the mode-mixing channel does not form part of the SU(1,1) group. Instead, such an interferometer would be described by a larger group, for example, the unitary group associated with a double covering of Sp($4,\mathbb{R}$).^{36} Similarly, our full interferometer with a general unitary Gaussian channel could be considered an Sp($6,\mathbb{R}$) interferometer. If we assume $N\xaf\u226b1$ in (8), we obtain $H\u2248A2\u2009sin2\u2009\theta (1\u2212sin2\u2009\theta \u2009sin2\u2009\varphi A)N\xafN/2$, where we have also assumed that $N\u226b1/2$ and taken $\nu A=\pi $.

### C. Sensitivity

We now choose a specific measurement process, the sum of the number of particles in the side modes. That is, the measured observable is $S\u0302=N\u0302:=a\u03021\u2020a\u03021+a\u03022\u2020a\u03022$. The square of the sensitivity of the interferometer is defined as (see, e.g., Ref. 37 for a derivation)

where $Var(S\u0302):=\u3008S\u03022\u3009\u2212\u3008S\u0302\u30092$. For Gaussian states, this is related to the Fisher information through^{38}

where $F0:=1/\Delta 2\epsilon $, such that $F\u2265F0$. Writing (9) in the CMF in the *q* and *p* basis, we find the simple expressions^{39}

where *n* is the number of modes, which is 2 in this case, and $\sigma s$ and $ds$ are the covariance and displacement matrices of the side modes, which are generated from the full ** d** and $\sigma $ by simply removing the first two rows and columns, respectively. Expressions for the mean $\u3008S\u0302\u3009$ and variance $Var(S\u0302)$ for the squeezing channel can be found in Appendix E when working with small

*s*. In the limit that $N\xaf\u226b1$ and $N\u226b1$, we find $F0=B2\u2009sin2(2\theta )NN\xaf/8$, where we have also taken the optimum phase relation $\nu B=0$. This matches the corresponding QFI, and so the number-sum measurement is an optimum measurement scheme in these limits.

When working with small *m*, full expressions for the mean $\u3008S\u0302\u3009$ and variance $Var(S\u0302)$ for the mode-mixing channel can be found in Appendix E. Taking $N\xaf\u226b1$ and $N\u226b1$, we obtain $F0\u2248A2\u2009sin2\u2009\theta (1\u2212sin2\u2009\theta \u2009sin2\u2009\varphi A)N\xafN/2$, where we have also taken the optimum phase relation $\nu A=\pi $. Therefore, as with the squeezing channel, this matches the QFI in the corresponding limits. Future work will apply multiparameter estimation techniques to, for example, estimate squeezing and mode-mixing parameters simultaneously.

## IV. IMPLEMENTATION IN A BEC

As an example of an implementation of the scheme, consider a one-dimensional BEC trapped in a box of length *L* and interacting with a periodic potential $\epsilon v\epsilon (x,t)$, with $\epsilon \u226a1$. The Hamiltonian for the Bose gas is given by

where $\Psi \u0302$ is the field operator of the atoms. The constant $g=4\pi \u210f2a/m$ is the coupling strength for a two-body contact potential, with *a* the s-wave scattering length and *m* the atomic mass.^{40,41} Depending on the nature of the potential $\epsilon v\epsilon (x,t)$, this implementation can be used to measure gravitational effects including gravitational fields and their gradient,^{42} gravitational waves,^{18} and also electromagnetic fields interacting with the phonons.^{43–45}

For a large number of atoms $N0\u226b1$ in the dilute and low-energy regime, the Bogoliubov approximation enables the decomposition of the atomic field operator $\Psi \u0302(x,t)=[\varphi 0a0+\psi \u0302(x,t)]e\u2212i\mu t/\u210f$ in terms of a classical function for the ground state $\varphi 0a0\u2248\varphi 0N0$ and a quantum field $\psi \u0302(x,t)$ that can be written as an infinite sum of quasiparticle modes, with *μ* the chemical potential. At low energies, these quasiparticles behave as phonons with frequency $\omega n=n\pi cs/L$, where *c _{s}* is the speed of sound. In Appendix H, we show how different choices of potentials $v\epsilon (x,t)$ and resonance conditions can implement $Ss(r),\u2009St(\theta )$, and $S\epsilon $ for the phonons.

Consider, for example, taking $v\epsilon (x,t)=v0x2\u2009sin\u2009\Omega t$, with $v0$ and Ω constants. If two phonon modes *l* and *n* have frequencies such that $\Omega =\omega n+\omega l$, then this potential will generate a resonant two-mode squeezing Gaussian channel that has squeezing parameter $s=\epsilon B4=2\epsilon |mlnvl,n|t/\u210f$, with

Here, $\zeta :=\u210f/(2mcs)$ is the BEC healing length, and the coefficients $vl,n$ are given by $vl,n=(1+(\u22121)l+n)v0$.^{42} If we assume that $Ss(r)$ creates $NP=2\u2009sinh2\u2009r$ phonons in modes *n* and *l*, and that the squeezing channel is generated by a quadratic potential, the precision in estimating *ε* is then

where we used $M=\tau /t$, with *τ* the integration time and *t* the interaction time, which is bounded by the lifetime of the phonons. Similarly, one can estimate the precision for the mode-mixing channel, where $A=8|mlnvl,n|t/\u210f$. However, in this case, *l* and *n* are mode numbers of frequency modes whose difference matches Ω.

Note that when the phonons interact with a periodic potential, which encodes the parameter to be estimated, it has been found to be convenient to choose phonon modes that are in resonance with the frequency of the potential. Previous work shows that resonances producing parametric amplification lead to higher degrees of mode-mixing and squeezing of phonons.^{29,31,42} The periodic potential can be due to gravitational^{1,42} or electromagnetic^{43–46} effects and lead to higher sensitivities in the detectors.^{17,18}

## V. APPLICATION OF FREQUENCY INTERFEROMETRY TO GRAVITATIONAL WAVE DETECTION

In the case that the external potential corresponds to a monochromatic continuous gravitational wave of frequency $\Omega =\omega l+\omega n$, the potential is $v\epsilon (x,t)=v0x2\u2009sin\u2009\Omega t$, with $v0=m\Omega 2/4$ (see Refs. 47 and 48 and Appendix J).

The precision (15) in the case of the estimation of the amplitude *ε* of the gravitational wave depends on multiple parameters that are inter-dependent. These include the phonon lifetime, the speed of sound, and the phonon mode numbers. It is convenient to write the precision in terms of the scattering length *a* and *α*, which is a constant that parameterizes the dimensions of the trap (see the Appendix I),

where *N _{d}* is the number of independent detectors. The parameters in the precision must satisfy constraints since we require the BEC to be a weakly interacting gas in the dilute regime ($N0a3/V\u226a1$, where

*V*is the volume of the BEC), and its excitations must be in the phononic regime ($\u210f\omega n,l\u226amcs2$). We work within the Bogoliubov approximation where the ground state operators can be treated classically and assumes that the ground state is macroscopically occupied in comparison to the excited states. We can approximate the number of excited atoms in a particular mode

*n*by $(mcs2/\u210f\omega n)Np$, and the total number of excited atoms must be significantly less than the number of condensed atoms for the Bogoliubov approximation to be satisfied. For a weakly interacting Bose gas, it then follows that the system can be described by a classical ground state with quantum quasiparticles in higher modes.

^{40}

It is out of the scope of this paper to provide optimal experimental parameters, which can only be determined using specialized numerical methods. The aim of this section is to illustrate the advantages of using frequency interferometry with the gravitational wave detector concept. For this reason, we only provide an example of “best case” parameters that give high sensitivities while satisfying all the constraints mentioned above. In Appendix I, we provide a table with a few extra examples and give further detail on the relationship between the different parameters and the constraints that they must satisfy. In Ref. 49, a comprehensive review of the gravitational wave detector will be provided with a detailed account of the optimal sensitivities that can be achieved using frequency interferometry, and with a range of constrained experimental parameters.

Our example considers $7Li$ BEC with $N0=4.4\xd7108$ atoms^{50} in a trap with $L=0.005\u2009m$. We assume that the scattering length is tuned to $a=99a0$, with *a*_{0} the Bohr radius. Taking into consideration two and three-body losses, we initially prepare the modes *n *=* *260 and *l *=* *258 in a two-mode squeezed state with *N _{p}* = 4000 phonons and a lifetime of $t=1\u2009ms$. We obtain a sensitivity of $\Delta \epsilon =4.3\xd710\u221221$ for the interferometry scheme at a gravitational wave of angular frequency $\Omega =25.3\u2009kHz$ and assuming ten independent detectors operating for $\tau =1\u2009yr$. In contrast, the previous non-interferometric scheme in Ref. 18 would have provided a sensitivity of $\Delta \epsilon =10\u221210$.

Note that we have chosen large mode numbers because they increase the precision; however, these numbers satisfy theoretically all the constraints. In practice, large mode numbers might be difficult to populate and to resolve (see discussion in Appendix I). We have chosen a phonon lifetime that is much shorter than the expected BEC and phononic lifetimes from three-body and two-body (Beliaev and Landau) decay processes for a three-dimensional BEC.^{51} This is because phonon decoherence would be expected to happen on scales shorter than the lifetimes. Quantum decoherence (or more strictly, dephasing) of phonons due to Beliaev and Landau damping as well as three-body atomic decay was considered in Refs. 51 and 52 within the Born–Markov approximation. It was found that very high squeezed states can decohere quickly due to Beliaev and three-body decay, potentially placing limitations on the original scheme.^{18} Here, we have used frequency interferometry to facilitate sensitivities with much smaller squeezing of phonons, and, thus, smaller decoherence rates. Also, three-body and two-body decay lifetimes and decoherence rates would be expected to very much set an upper bound for the BEC of the detector since it is assumed to operate in the quasi-one-dimensional or pure one-dimensional regimes, where decay processes are heavily suppressed.^{51,52} Therefore, a more comprehensive review of decoherence is needed, concentrating on one-dimensional BECs, which will be performed in future work.^{49} The effect of temperature on the original scheme of the gravitational wave detector was also analyzed in Ref. 19. It was found that this had a negligible effect on the sensitivity of the detector. In addition to the effects of temperature and electromagnetic self-interactions in the BEC, other sources of noise for the detector will include Newtonian noise as is the case with LIGO, where a precise suspension mechanism is used to strongly suppress this, but which will be low in the very high frequencies searched for by this detector, and noise from fluctuations in the trapping potential of the BEC, which can be made very low with the use of permanent magnets.

The aim of this example is to show that higher sensitivities can be reached because frequency interferometry increases the number of phonons in the side modes by beam-splitting the phonons with the condensate, which can be achieved by, for example, modifying the trapping potential (see Refs. 53 and 54 and Appendix H). For a more detailed derivation, constraints, and examples on sensitivities, see Appendix I. In Appendix H 2, we outline a method for creating squeezed states of phonons.^{45} This method predicts theoretically that for the parameters considered in the example, *N _{p}* = 4000 phonons can be squeezed, which corresponds to a squeezing strength of

*r*=

*4.5. Further investigation is required to confirm whether this could be obtained in practice. A squeezed state of phonons was attempted in the laboratory*

^{46}using a similar method; however, it was not possible to determine if the state had quantum properties. A two-mode quantum squeezed state of phonons with a low level of squeezing (reaching a squeezing strength of around,

^{55}

*r*=

*0.9) was claimed to have been created in an experiment looking for Hawking radiation of analogue black holes and its entanglement.*

^{56}

The theoretical proposal of using BEC phonons to detect high-frequency gravitational waves was introduced in Ref. 18 considering one-dimensional traps. The effect of gravitational waves on phonons in three-dimensional BECs has been studied for free falling^{57,58} and rigid boundary conditions.^{59} The authors in Refs. 52, 57, and 60 expressed concerns about the high level of squeezing and long phonon lifetimes required to detect gravitational waves in realistic experiments. In this paper, we show that these requirements can be dramatically reduced using frequency interferometry where the condensate is mixed with the phonon modes using a tritter. Mixing the condensate with phonon modes was suggested in Ref. 48. Here, we implement this idea concretely using frequency interferometry and prove mathematically that the sensitivity is, indeed, improved. Other ideas have been suggested, including using inhomogeneous condensate flows, such as a vortex lattice.^{60} In Ref. 61, the authors argued that modulating the trap potential induces parametric resonances in the BEC, which could lead to an enhancement in detection. Parametric amplification is known to be equivalent to squeezing in the quantum regime, see, e.g., Ref. 42 for this connection in BECs where an oscillating gravitational field is shown to create parametric amplification and is equivalent to squeezing the phonons of the BEC. Gravitational wave detection in cold atoms that exploits the interference of spacetime separated atomic matter waves, as opposed to the interference of phonon modes in a strongly interacting BEC, has also been studied, for example, in Refs. 62–65. Furthermore, superfluids have also been considered for gravitational wave detection.^{66,67}

There are no known astrophysical objects that are small and dense enough to emit at frequencies beyond $104\u2009Hz$.^{68} However, we consider sensitivities to higher frequencies because any discovery beyond this range would be produced by other yet-unknown sources of gravitational waves in the cosmos. This includes exotic objects, such as primordial black holes or boson stars, and cosmological events in the early universe, such as phase transitions, preheating after inflation, oscillons, and cosmic strings.^{68,69} Detecting gravitational waves at frequencies beyond the range of sensitivity of LIGO might eventually lead to the detection of dark matter. Interestingly, ultralight dark matter models for collective phenomena, as opposed to single scattering regimes, predict monochromatic long-lived dark matter waves at very high frequencies. The frequency of these wave is set by the dark matter mass and ranges from $10\u22128$ to 10^{14} Hz.^{70} The range of the application of frequency interferometry to the detection of long-lived gravitational waves depends on the frequency of the modes that resonate with the wave. In the case that the modes are BEC phonons, this is given by the mode number, the speed of sound, and the dimensions of the trap.

Interestingly, Weber bars attempted to use resonances of frequency modes to detect gravitational waves. However, the large metallic devices were not cooled below $100\u2009mK$. At these temperatures, the vibrational modes are in the classical domain and cannot be prepared in quantum states. Ultracryogenic detectors could reach the standard quantum limit, but squeezing of the modes could not be produced in these systems.^{71} The proposal in Ref. 18 resembles a quantum version of a Weber bar, where the atom–atom interactions in the BEC are Hamiltonian non-linearities, which produce quantum excitations of phononic modes. In this case, a harmonic perturbation can produce squeezing via parametric amplification exploiting resonances between the potential perturbation and the quantum modes.

## VI. COMPATIBILITY WITH GENERAL RELATIVITY

Quantum spatial atom interferometers are commonly described by non-relativistic quantum mechanics. The state evolution is then given by the Schrödinger equation, which is invariant under Galilean transformations. For this reason, the description is only compatible with the Newtonian approximation of gravity where the notion of time is absolute. Describing spatial interferometers beyond the Newtonian approximation is non-trivial since the equations as well as the inner products must be Lorentz-invariant and conserve quantum probabilities. This condition can be consistently satisfied for quantum fields but is problematic for individual particles.^{72} While a few studies have made progress toward describing the quantum states of single particles in curved spacetime,^{73–75} a consistent covariant description of atom interferometry still remains an open question.

Phonons are excitations of a massless quantum field, where the speed of propagation is the speed of sound. For this reason, it is natural to use the mathematical formalism of quantum field theory in curved spacetime in the frequency interferometry proposal introduced here. Frequency interferometry can then be applied in both Newtonian and General Relativistic regimes, in the latter using a theoretical description that is underpinned by quantum field theory in curved spacetime.^{44,76} This enables the application of the scheme to study effects in General Relativity and modified theories, including the estimation of spacetime parameters, and searches for dark energy and matter.

## VII. SUMMARY

We have introduced a scheme for interferometry that uses the frequency modes of a quantum field trapped in a localized potential. This setup does not require a large spatial extent, facilitating miniaturization. Since the modes occupy the same region in space at all times, it is possible to estimate parameters encoded in both phase-shift and global channels acting on the modes. This includes channels that entangle the modes. As an example, we estimate the parameters of two-mode Gaussian channels using a three-mode scheme that implements an analogue in the frequency domain of the tritter operation introduced in the pumped-up SU(1,1) scheme.^{22} The tritter operation improves the scaling of the precision with the number of particles. Quantum frequency interferometry should be implementable in many systems, such as optical, hybrid atom-light, superconducting circuits, cold atoms, and BECs, where the global channels could be generated using, for example, non-linear mediums and non-linear interactions. We show how the three-mode example can be implemented in a BEC, where the scheme exploits atom–atom interactions. These type of interactions produce undesired noise in other setups such as atom spatial interferometry. When the external potential that interacts with the BEC corresponds to an electromagnetic or a gravitational field, the scheme can be used to estimate the field parameters with a system trapped in a millimeter-scale trap. We show that this scheme can improve the precision in the detection of gravitational waves, enabling good sensitivities even in the case that squeezing is much smaller than assumed previously in Ref. 18, and that the system suffers from short phononic lifetimes.

## DEDICATION

We would like to dedicate this paper to Roger Penrose in celebration of his Nobel Prize and 90th Birthday. We found great inspiration in his work, and we are grateful for everything we have learned from him during the many discussions we were fortunate to share with him.

## ACKNOWLEDGMENTS

We thank Paul Juschitz, Jan Kohlrus, Daniel Goldwater, Tupac Bravo, Daniel Hartley, and Dennis Rätzel for useful discussions and comments. R.H. and I.F. would like to acknowledge that this project was made possible through the support of a donation by John Moussouris and the grant “Leaps in cosmology: gravitational wave detection with quantum systems” (No. 58745) from the John Templeton Foundation. R.H. would also like to acknowledge the support of the ID 61466 grant from the John Templeton Foundation, as part of the QISS project. The opinions expressed in this publication are those of the authors and do not necessarily reflect the views of the John Templeton Foundation.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Richard Howl:** Conceptualization (supporting); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Ivette Fuentes:** Conceptualization (lead); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (lead); Resources (lead); Supervision (lead); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

### APPENDIX A: SYMPLECTIC MATRICES OF INTERFEROMETRY OPERATIONS

Here, we provide the symplectic matrices, in the real *q* and *p* representation, for the various processes involved in our considered active interferometry schemes. As in the main text, the initial state of the pump mode is assumed to be a coherent state, and so the displacement and covariance matrices of the full input state to the interferometer are

where we have written $\alpha \u2261N\xafei\u03d10$, with $N\xaf$ the total particle number and $16:=diag(1,1,1,1,1,1)$ the identity matrix, of which the first two rows and columns are for the pump mode, the next two rows and columns are for one of the side modes, and the final rows and columns are for the other side mode.

The first stage of the interferometer is the two-mode squeezing operation that parametrically populates the side modes, which has the following symplectic matrix (see, e.g., Ref. 25),

where *r* is the squeezing parameter and $\u03d1sq$ is the squeezing phase. Here, the first two columns and rows are for the pump, the next two column and rows are for one of the side modes, and the last two columns and rows are for the other side mode.

The next stage is a tritter between the side modes and the pump, which has the symplectic matrix given in (B11). Following the tritter, there is the squeezing or mode-mixing channel, which are defined by the unitary transformations:

or

where $\xi :=sei\varphi B$ and $\zeta :=mei\varphi A$, with $s\u22650,\u2009m\u22650$, and $\varphi A,\varphi B\u2208\mathbb{R}$. The symplectic matrices for these unitary evolutions are (see, e.g., Ref. 25)

and

where

In contrast, the symplectic matrix for a unitary phase evolution $U\u0302(\varphi )=exp\u2009(\u2212i\varphi N\u0302/2)$ would be the following (see, e.g., Ref. 25):

### APPENDIX B: DERIVATION OF SYMPLECTIC MATRIX OF TRITTER

The tritter used in the pumped-up SU(1,1) interferometry scheme is generated by the following Hamiltonian:^{23}

which, in the Heisenberg picture, results in

where $\theta :=Gt$ is the angle for an evolution time *t* and ϑ is the phase. We can write the above transformation in the matrix form

where

We now move to the real symplectic *q* and *p* representation,

where

and

The symplectic representation of the tritter transformation is then

where

Note that with a conventional two-way beam splitter, $\theta =\pi /2$ would swap the modes. However, for the above tritter, $\theta =\pi /2$ would not completely swap the side modes and pump modes. This is responsible for *N* appearing in the QFI expressions even when $\theta =\pi /2$.

### APPENDIX C: QUANTUM FISHER INFORMATION OF PHASE-SHIFT CHANNEL

If a conventional phase-shift channel were used instead of the squeezing or mode-mixing channels, the QFI is

where

### APPENDIX D: TURNING POINTS OF QUANTUM FISHER INFORMATION

As discussed in the main text, the QFI *H* of the interferometer with a squeezing channel has three turning points when the optimum phase relations are chosen: $\u03d1sq=\varphi B+\pi /2$ and $\u03d1=\u03d10\u2212\varphi B/2+\pi /4$. These occur at *θ* = 0, $\theta =\pi /2$, and $\theta =\theta t$, where $\theta t=cos\u22121(zt)/2$ with

The angle *θ _{t}* matches that found in Ref. 23 for a phase-shift channel. When $N\xaf$ is large, it can be approximated by

^{23}

For θ = 0, $\theta =\pi /2$, and $\theta =\theta t$, H is

where, for the $\theta =\theta t$ turning point, we have assumed that $N\xaf\u226b1$ and taken the optimum phase relation $\u03d1=\u03d10+\u03d1sq/2\u22122\varphi B$.

Analogous expressions to (D4)–(D6) for the squeezing case can be obtained for *H* of the mode-mixing channel when *θ* = 0 and $\theta =\pi /2$, with $\varphi A=\pi /2$, and $\theta =\pi /2$ with $\varphi A=0$,

where, for the last case, we have assumed that $N\xaf\u226b1$ and taken $\u03d1=\u03d10\u2212\u03d1sq/2+\pi /2$.

### APPENDIX E: SENSITIVITY

In the main text, we assume an intensity measurement. For the squeezing channel and working with small *s*, we find the following expressions for this measurement scheme:

where

We now use particle conservation to write $|\alpha 0(r)|2=|\alpha |2\u22122\u2009sinh2\u2009r$ and $|\alpha |2=N\xaf$, so that $|\alpha (r)|2=1/\epsilon \u22122\u2009sinh2\u2009r$, where $\epsilon :=1/N\xaf$, and take $N\xaf\u226b1$ to find

where in the last line, we have assumed that $N\u226b1$ and taken the optimum phase relation $\nu B=0$.

If, on the other hand, the mode-mixing channel is chosen, then, for small *m*, we have

where

Taking $N\xaf\u226b1$, we find

where in the last line, we have assumed that $N\u226b1$ and taken the optimum phase relation $\nu A=\pi $.

### APPENDIX F: FULL UNDEPLETED PUMP REGIME

In the main text, we have assumed that the pump is relatively undepleted after the first active element (see also Ref. 23). If we want to further assume that the pump is also relatively undepleted after the tritter stage, then *θ* cannot be too large. After the tritter stage, in general, the number of particles in the pump and side modes is the following:

Let us require that $N=\gamma N0$ and $N(\theta )=\delta N0(\theta )$, where $\gamma \u226a1,\u2009\delta \u226a1$, and $\delta \u2265\alpha $. Then, *θ* must satisfy

For example, taking $\delta =0.1$, we obtain

which, in the limit $\gamma \u21920$, gives $\theta 2\u22480.0938$, and we note that $\theta 2/\u2009sin2\u2009\theta \u22481.03$.

In the case of the squeezing channel, the QFI in this fully undepleted pump regime becomes

where we have again used $\u03d1=\u03d10+\u03d1sq/2\u22122\varphi B$ and further assumed that $r\u226b1$ and $N\xaf\u226b1$ in the last line. When considering a number-sum measurement scheme, the square-inverse of the sensitivity is

when $r\u226b1$, and $2\u03d1=2\u03d10+\u03d1sq\u22122\varphi B$. This agrees with the QFI expression above (F6). The number-sum measurement is, therefore, an optimum measurement scheme in these limits.

For the pump to remain relatively undepleted before the mode-mixing channel, the QFI becomes

where we have assumed that $r\u226b1$ in the last line. This is similar to the QFI for the squeezing channel (F5), just with *B* replaced by *A*. Equally, the square-inverse of the sensitivity is

with $2\u03d1=2\u03d10\u2212\u03d1sq$ and $r\u226b1$.

### APPENDIX G: HETERODYNE DETECTION

Rather than using a number-sum measurement, another possibility would be to use a heterodyne measurement, for example, between the pump and the side modes. Balanced homodyne detection for the side modes was considered in Ref. 78 for a standard SU(1,1) interferometer and in Ref. 79 for a “truncated” SU(1,1) experiment. In our considered heterodyne case, at the measurement stage, a balanced beam splitter could be applied between one of the side modes and the pump, and the difference of the number of particles in the two output parts of the final beam splitter could be considered: $S\u0302=N\u03021\u2212N\u03022$. In the covariance matrix formalism, we have

where

However, in order to measure the squeezing parameter of the estimation channel, $Var(S\u0302)$ or $\u3008S\u03022\u3009$ would need to be considered as the signal (see, e.g., Ref. 80), and, therefore, the variance of this would be used in the error estimation.

### APPENDIX H: IMPLEMENTATION IN A BEC

In this section, we provide detail on the implementation of the three-mode frequency interferometric scheme using phonons of a BEC. The Hamiltonian of a Bose gas with a trapping potential $v(r)$ is

where *V* is the volume. We consider the *dilute regime*, where the contact potential only depends on two-body interactions and the coupling strength is given by $g=4\pi \u210f2a/m$, where *a* is the s-wave scattering length.^{40,41} It is convenient to decompose the field operator $\Psi \u0302\u2020$ in terms of the annihilation operator for the ground state $a\u03020$ and the operators for the *n*th excited state $a\u0302n$,

where

with the commutator $[a\u0302n,a\u0302l\u2020]=\delta nl$. We consider $\mu (t)=\mu t$, where *μ* is the chemical potential. In the case that the external potential vanishes, the Hamiltonian can be diagonalized by taking two steps. The first is assuming that the ground state is macroscopically occupied in a large coherent state. In this case, the ground state operator can be replaced by $a\u03020\u2248N0$. This is called the *Bogoliubov approximation*, and it is applicable to systems with a large number of particles in the low temperature regime (*T* much smaller than the condensates' critical temperature). The second step involves applying a Bogoliubov transformation to $a\u0302n$, such that

where $[b\u0302n,b\u0302l\u2020]=\delta nl$. Neglecting trilinear and quartic terms in $b\u0302n,b\u0302n\u2020$ (since they have fewer factors of $N0\u226b1$) yields the diagonal Hamiltonian (see, e.g., Ref. 40),

Here $::$ refers to normal ordering and *u _{n}* and

*v*are mode solutions to the Bogoliubov–de-Gennes equations,

_{n}^{40}

satisfying the orthonormal condition

The ground state wave function $\varphi 0$ satisfies the time-independent Gross–Pitaevskii equation

such that $\varphi 0(t):=\varphi 0e\u2212i\mu t/\u210f$ satisfies the time-dependent version. The energy spectrum is given by

where $cs:=g\rho /m$ is the speed of sound of the BEC and *ρ* is the number density. In the *low-energy limit*$\u210f\omega n\u226amcs2$, the dispersion law is linear and, thus, $b\u0302n\u2020$ and $b\u0302n$ create and annihilate phonons of the BEC.

We now apply a small time-dependent potential $\epsilon v\epsilon (r,t)$ to the BEC where $\epsilon \u226a1$. This introduces a term $\epsilon v\epsilon \Psi \u0302\u2020\Psi \u0302$ to (H1), which, after applying (H2) and (H5), provides an interaction Hamiltonian (see the Appendix C of Ref. 42 for a detailed derivation using the grand canonical Hamiltonian),

where $exp\u2009(i\u03d1(r)):=\varphi 0(r)un\u2217(r)+\varphi 0\u2217(r)vn\u2217(r)$.

We now consider how the experimentalist can tailor a sequence of external potential of the form $v\epsilon (r,t)$ with the appropriate resonance conditions to implement the unitaries $U\u0302(\theta )$ and $U\u0302(-\theta )$ in the phases (i) and (iii) of our frequency interferometric scheme. Frequency interferometry can then be performed with a single BEC in a pumped-up SU(1,1) scheme. In this case, the condensate atoms act as the pump, and two phonon modes can be used as the side modes.

##### 1. Tritter

The Hamiltonian for a tritter is given by (B1) where here we treat $a\u03020$ as the annihilation operator for the condensate, and $a\u0302n\u22600$ as the annihilation operator for the phonon modes, which we denoted as $b\u0302n\u22600$ above. Since the condensate must be more populated than the phonon modes before and after the tritter for our description of the BEC used above to still hold, we can apply the Bogoliubov approximation and drop the hat on $a\u03020$, leaving us with

This can be picked out from (H12) by choosing an oscillating potential of the form $V(t)=\epsilon V0\u2009cos\u2009(\Omega t)\u2009cos\u2009(\Omega \u2032t)$, where $\Omega :=\omega m+\omega n$ and $\Omega \u2032:=\omega n\u2212\omega n$, and assuming that $\u03d1n(r)\u2248\u03d1m(r)$, which could be achieved, for example, by choosing modes with equal and opposite momenta in a uniform BEC with periodic boundary conditions.^{40}

##### 2. Two-mode squeezing

To create a two-mode squeezed state of phonons, we require a Hamiltonian of the form (see, e.g., Ref. 25)

This can be obtained from (H12) by choosing an oscillating potential to pick out these particular terms on resonance.^{42} For example, $V\epsilon (t)=\epsilon V0\u2009sin\u2009\Omega t$ would achieve this, where $\Omega :=\omega m+\omega n$ and *V*_{0} is a constant amplitude. The amount of squeezing that can be generated for phonons by periodically shaking the trap as proposed here has been considered for existing experimental setups, finding that high levels of squeezing are theoretically possible,^{45} even when taking into account decoherence due to phonon interactions.^{51}

### APPENDIX I: APPLICATION TO THE DETECTION OF GRAVITATIONAL WAVES

As in the previous appendix, we assume a Bose gas operating in the dilute regime. However, we now take a box potential with $v(r)=0$ and $V=La$, where *L* is the length and $a$ is the cross-section. We assume that *L* is much greater than the dimensions of $a$ and restrict the analysis to modes with vanishing transversal wave numbers, i.e., we only consider the direction along *L*, which we label *x*. The mode solutions, which fulfill von Neumann boundary conditions since *ρ* vanishes at the potential walls, are given by^{40,42}

where

and $kn=n\pi /L=\omega n/cs$. Taking the time-dependent potential $v\epsilon (r,t)=v\epsilon (x,t)=v\epsilon (x)\u2009sin\u2009\Omega t$, writing $sin\u2009\Omega t\u22611/(2i)(ei\Omega t\u2212e\u2212i\Omega t)$, and applying the rotating-wave approximation, we yield^{42}

with

and $\zeta :=\u210f/(2mcs)$ is the healing length. If the external potential is linear, $v\epsilon (x)=v0x$, where $v0$ is a constant, then the coefficients are given by $vl,n=\u2212(1\u2212(\u22121)l+n)(v0/L)$, and if it is quadratic $v\epsilon (x)=v0x2$, then $vl,n=(1+(\u22121)l+n)v0$.^{42}

From inspection of the interaction Hamiltonian, we see that the resonance condition $\Omega =\omega n$ generates a displacement of the mode *n*, which creates a (classical) coherent state. Two-mode squeezing (parametric amplification) is produced through the resonance condition $\Omega =\omega n+\omega m$ and $\Omega =\omega n\u2212\omega m$ leads to mode-mixing (frequency conversion). In the two-mode squeezing case, the squeezing parameter is $s=\epsilon B4=2\epsilon |mlnvl,n|t/\u210f$ where

with *l* and *n* ($l\u2260n$) the mode numbers of the two frequency modes whose sum is resonant with Ω. For the mode-mixing channel, $A=8|mlnvl,n|t/\u210f$; however, *l* and *n* are mode numbers of frequency modes whose difference sum is resonant with Ω.

Using the results also from the previous appendix, we can perform frequency interferometry with a single BEC in a pumped-up SU(1,1) scheme, with the condensate atom acting as the pump and two phonon modes as the side modes. Since the number of condensate atoms *N*_{0} must be much larger than the number of phonons *N _{p}*, with the implementation of a number-sum measurement scheme, the sensitivity of the interferometer to a two-mode squeezing channel is

where *M* is the number of repetitions of estimation and *θ* must satisfy (F3). In the second line, we have assumed that $\theta 2N0Np\u226bNP2$, and in the last line, we have assumed that $M=\tau /t$, where *τ* is the total time of the full estimation procedure, with *t* then being the time it takes for each individual estimation (which we take to be the phonon lifetime). *N _{d}* is the number of independent detectors.

We now consider a gravitational wave interacting with a BEC. Gravitational waves are often considered from the perspective of the transverse-traceless (TT) frame. However, this is not the frame of a BEC experimentalist, which can lead to intuitive effects that must be considered carefully, such as the fact that a rigid cavity oscillates in this frame. Instead, we consider the BEC from the proper detector frame,^{47,48} which is the frame closest to a BEC experimentalist (for example, rigid cavities do not appear to change length in this frame). In this frame, a GW just acts as a Newtonian force obeying Hooke's law:^{47} $Fi=mh\xa8ijTT\xi j/2$, where *ξ ^{j}* is the coordinate separation,

*m*is the mass, and $hijTT$ is the spatial gravitational wave perturbation in the

*TT*frame. For a monochromatic wave of frequency Ω, $hijTT=\epsilon \u2009sin\u2009(\Sigma t)$

^{47}and so $Fi=m\Omega 2\u2009sin\u2009(\Sigma t)\xi j/2$. This results in a quadratic potential term on a quasi-one-dimensional BEC with $v\epsilon (x,t)=mx2\Omega 2\u2009sin\u2009\Omega t/4$. A full detailed derivation of this potential starting from the Lagrangian of the BEC is given in Appendix J. Comparing to the above, the gravitational wave will create parametric-amplification in the BEC, leading to squeezing of phonons.

^{48}Using (I12) and taking $\Omega =\omega l+\omega n$ with

*l*+

*n*even, we find

where we have used $cs=g\rho /m$ with $\rho =N0/V$, defined $a:=\pi R2$ with $R=:\alpha L$, and *N _{d}* is the number of independent detectors. This improves on sensitivity of Ref. 18, which, written in fundamental experimental parameters, would be

Note that the sensitivity reported in Ref. 18 is given in terms of $\Delta \epsilon /\Omega $. However, we use $\Delta \epsilon $ since it is a figure of merit, which is more convenient for resonant detectors. The sensitivity reported in Ref. 18 also has a different dependence on *N*_{0}. This is because a factor $1/N0$ was added, assuming that the sensitivity scales with the number of atoms detected in the experiment. In this work, the dependence on *N*_{0} is not assumed, it is obtained from first principles calculations.

Taking, for example, an $7Li$ BEC with $N0=4.4\xd7108$, *N _{p}* = 4000, $a=99a0$ (

*a*

_{0}is the Bohr radius), $L=0.005\u2009m,\u2009\alpha =0.001$, and $t=1\u2009ms$, and resonating with modes

*n*=

*260 and*

*l*=

*258, we obtain a sensitivity of $\Delta \epsilon =4.3\xd710\u221221$ for the interferometry scheme at a gravitational wave of angular frequency $\Omega =25.3\u2009kHz$ and assuming ten independent detectors operating for $\tau =1\u2009yr$. In contrast, the previous non-interferometric scheme would have provided a sensitivity of $\Delta \epsilon =10\u221210$. In Tables I and II, we show further examples of the precision reached by the phononic gravitational wave detector using $7Li$ (in the $|F=1,mF=1\u3009$ hyperfine state) and $87Rb$ BECs. The frequency interferometry (FI) scheme improves by many orders of magnitude of the sensitivities, given by Ref. 18. The integration time*

*τ*is one year, the tritter angle is $\theta =0.31$, and we considered

*N*= 10 independent detectors. For $7Li$, the three-body loss rate is calculated using the universal theory of Efimov physics,

_{d}^{81,82}which have been shown to match experimental data very well.

^{83}For $87Rb$, we assume a three-body rate coefficient of $4\xd710\u221230\u2009cm6\u2009s\u22121$.

^{84}

Parameter Units . | L mm
. | α
. | N_{0}
. | N
. _{p} | l
. | n
. | t s
. | $a/a0$ . | n_{0} $cm\u22123$
. | Ω kHz . | c m/s
. _{s} | $\Delta \epsilon $ . | $\Delta \epsilon /\Omega $ $Hz\u22121/2$ . |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|

Original scheme | 4 | 0.01 | $4.4\xd7109$ | 2500 | 480 | 478 | 0.02 | 119 | $2.2\xd71014$ | 28.4 | 0.04 | $4.1\xd710\u221211$ | $2.4\xd710\u221213$ |

FI scheme | 4 | 0.01 | $4.4\xd7109$ | 2500 | 480 | 478 | 0.02 | 119 | $2.2\xd71014$ | 28.4 | 0.04 | $2.9\xd710\u221221$ | |

Original scheme | 5 | 0.001 | $4.4\xd7108$ | 4000 | 260 | 258 | 0.001 | 99 | $1.1\xd71015$ | 25.3 | 0.08 | $1.0\xd710\u221210$ | $6.4\xd710\u221213$ |

FI scheme | 5 | 0.001 | $4.4\xd7108$ | 4000 | 260 | 258 | 0.001 | 99 | $1.1\xd71015$ | 25.3 | 0.08 | $4.3\xd710\u221221$ | |

Original scheme | 4 | 0.01 | $4.4\xd7109$ | $4.4\xd7107$ | 480 | 478 | 0.04 | 99 | $2.2\xd71014$ | 25.9 | 0.03 | $1.8\xd710\u221215$ | $1.1\xd710\u221217$ |

FI scheme | 4 | 0.01 | $4.4\xd7109$ | $4.4\xd7107$ | 480 | 478 | 0.04 | 99 | $2.2\xd71014$ | 25.9 | 0.03 | $2.0\xd710\u221223$ |

Parameter Units . | L mm
. | α
. | N_{0}
. | N
. _{p} | l
. | n
. | t s
. | $a/a0$ . | n_{0} $cm\u22123$
. | Ω kHz . | c m/s
. _{s} | $\Delta \epsilon $ . | $\Delta \epsilon /\Omega $ $Hz\u22121/2$ . |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|

Original scheme | 4 | 0.01 | $4.4\xd7109$ | 2500 | 480 | 478 | 0.02 | 119 | $2.2\xd71014$ | 28.4 | 0.04 | $4.1\xd710\u221211$ | $2.4\xd710\u221213$ |

FI scheme | 4 | 0.01 | $4.4\xd7109$ | 2500 | 480 | 478 | 0.02 | 119 | $2.2\xd71014$ | 28.4 | 0.04 | $2.9\xd710\u221221$ | |

Original scheme | 5 | 0.001 | $4.4\xd7108$ | 4000 | 260 | 258 | 0.001 | 99 | $1.1\xd71015$ | 25.3 | 0.08 | $1.0\xd710\u221210$ | $6.4\xd710\u221213$ |

FI scheme | 5 | 0.001 | $4.4\xd7108$ | 4000 | 260 | 258 | 0.001 | 99 | $1.1\xd71015$ | 25.3 | 0.08 | $4.3\xd710\u221221$ | |

Original scheme | 4 | 0.01 | $4.4\xd7109$ | $4.4\xd7107$ | 480 | 478 | 0.04 | 99 | $2.2\xd71014$ | 25.9 | 0.03 | $1.8\xd710\u221215$ | $1.1\xd710\u221217$ |

FI scheme | 4 | 0.01 | $4.4\xd7109$ | $4.4\xd7107$ | 480 | 478 | 0.04 | 99 | $2.2\xd71014$ | 25.9 | 0.03 | $2.0\xd710\u221223$ |

Parameter Units . | L mm
. | α
. | N_{0}
. | N
. _{p} | l
. | n
. | t s
. | $a/a0$ . | n_{0} cm^{−3}
. | Ω kHz . | c m/s
. _{s} | $\Delta \epsilon $ . | $\Delta \epsilon /\Omega $ $Hz\u22121/2$ . |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|

Original scheme | 6 | 0.0005 | $4.4\xd7108$ | 4400 | 480 | 478 | 0.002 | 109.6 | $2.6\xd71015$ | 5.3 | 0.01 | $3.8\xd710\u221210$ | $5.3\xd710\u221212$ |

FI scheme | 6 | 0.0005 | $4.4\xd7108$ | 4400 | 480 | 478 | 0.002 | 109.6 | $2.6\xd71015$ | 5.3 | 0.01 | $4.5\xd710\u221221$ | |

Original scheme | 4 | 0.003 | $4.4\xd7109$ | $4.4\xd7107$ | 1680 | 1678 | 0.003 | 109.6 | $2.4\xd71015$ | 25.65 | 0.01 | $7\xd710\u221215$ | $4.4\xd710\u221217$ |

FI scheme | 4 | 0.003 | $4.4\xd7109$ | $4.4\xd7107$ | 1680 | 1678 | 0.003 | 109.6 | $2.4\xd71015$ | 25.65 | 0.01 | $6.3\xd710\u221224$ |

Parameter Units . | L mm
. | α
. | N_{0}
. | N
. _{p} | l
. | n
. | t s
. | $a/a0$ . | n_{0} cm^{−3}
. | Ω kHz . | c m/s
. _{s} | $\Delta \epsilon $ . | $\Delta \epsilon /\Omega $ $Hz\u22121/2$ . |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|

Original scheme | 6 | 0.0005 | $4.4\xd7108$ | 4400 | 480 | 478 | 0.002 | 109.6 | $2.6\xd71015$ | 5.3 | 0.01 | $3.8\xd710\u221210$ | $5.3\xd710\u221212$ |

FI scheme | 6 | 0.0005 | $4.4\xd7108$ | 4400 | 480 | 478 | 0.002 | 109.6 | $2.6\xd71015$ | 5.3 | 0.01 | $4.5\xd710\u221221$ | |

Original scheme | 4 | 0.003 | $4.4\xd7109$ | $4.4\xd7107$ | 1680 | 1678 | 0.003 | 109.6 | $2.4\xd71015$ | 25.65 | 0.01 | $7\xd710\u221215$ | $4.4\xd710\u221217$ |

FI scheme | 4 | 0.003 | $4.4\xd7109$ | $4.4\xd7107$ | 1680 | 1678 | 0.003 | 109.6 | $2.4\xd71015$ | 25.65 | 0.01 | $6.3\xd710\u221224$ |

In deriving these results, we have been careful to satisfy various theoretical and experimental constraints. For example, we have made sure that the chosen modes are phononic modes ($\u210f\omega n,l\u226amcs2$), that we are operating in the dilute regime ($N0a3/V\u226a1$, where *V* is the volume of the BEC), that (F3) is satisfied, and that the Bogoliubov approximation is satisfied. The latter approximation is that the ground state operators can be treated classically and assumes that the ground state is macroscopically occupied in comparison to the excited states. For a weakly interacting Bose gas, it then follows that the system can be described by a classical ground state with quantum quasiparticles in higher modes.^{40} For modes *n* where $\u210f\omega n\u226amcs2$, the quasiparticles behave as phonons, which consist of multiple excited atoms. The average number of atoms making up a phonon in a particular state can be obtained from the Bogoliubov transformation that relates the atomic and quasiparticle bases.^{40} We can then approximate the number of excited atoms in a particular mode *n* by $(mcs2/\u210f\omega n)Np$, and the total number of excited atoms must be significantly less than the number of condensed atoms for the Bogoliubov approximation to be satisfied. In the cases of very high phonon numbers $Np\u22484.4\xd7107$ in Tables I and II, the number of excited atoms is about 0.1% of the number of condensate atoms. Here, although theoretically the Bogoliubov approximation is still satisfied, it is experimentally possible that effects beyond the Bogoliubov approximation could be seen.

In Tables I and II, we have assumed very low temperatures, so that there is macroscopic condensation. At very low temperatures, Landau damping of phonons can be neglected, which depends strongly on temperature.^{51} However, there is still two-body Beliaev decay as well as three-body atomic decay, which are approximately independent of temperature. Quantum decoherence (or more strictly, dephasing) of phonons due to Beliaev and Landau damping as well as three-body atomic decay was considered in Refs. 51 and 52 in the Born–Markov approximation. It was found that very high squeezed states can suffer a loss of squeezing very quickly due to Beliaev and three-body decay, potentially placing limitations on the original scheme.^{18} Here, we have used frequency interferometry to facilitate sensitivities with much smaller squeezing of phonons and, thus, smaller decoherence rates. We have also chosen phonon lifetimes that are much smaller than the expected ones from Beliaev and three-body decay in order to account for the fact that decoherence of squeezed states will occur faster than the predicted damping lifetimes. However, a more comprehensive review of decoherence for the assumed measurement schemes is needed, which will be performed in future work.^{49} This will take into account that loss of phonon squeezing does not occur in an exponential fashion—see Ref. 51, where, although there can be a sudden rapid loss of squeezing, large values of squeezing can still persist long after this sudden loss and thus effective decoherence time scales. Therefore, using decoherence times in estimating the sensitivity can be misleading. It should also be noted that three-body and two-body decay lifetimes and decoherence rates would be expected to very much set an upper bound for the BEC of the detector since it is assumed to operate in the quasi one-dimensional or pure one-dimensional regimes where decay processes are heavily suppressed.^{51,52} Furthermore, the effect of temperature on the original scheme of the gravitational wave detector was also analyzed in Ref. 19, where it was considered how a non-zero temperature will result in a partially mixed rather than pure state of phonons. It was found that this had a negligible effect on the sensitivity of the detector.

In addition to the effects of temperature and electromagnetic self-interactions in the BEC, other sources of noise for the detector will include Newtonian noise as is the case with LIGO, where a precise suspension mechanism is used to strongly suppress this, but which will be low in the very high frequencies searched for by this detector, and noise from fluctuations in the trapping potential of the BEC, which can be made very low with the use of permanent magnets. A thorough analysis of noise in the detector will be performed in Ref. 49.

The aim of Tables I and II is to provide a sense of what are the physical requirements and sensitivities involved. We have chosen parameters that satisfy all the constraints described above and that are, in principle, accessible to current experiments. We choose best case phonon numbers by considering that high numbers but making sure that the numbers chosen satisfy all the constraints mentioned above. Choosing high phonon numbers improves sensitivities; however, the use of such modes needs further consideration. While the scheme relies only on the resolution of squeezing, the frequency modes themselves cannot be resolved in a single direct measurement on the phonons due to the small phonon lifetimes. Resolving the modes requires integrating over many measurements^{46} or using indirect measurement schemes that prolong the frequency modes beyond the phonon lifetimes, such as phonon evaporation.^{53} While we conjecture that this does not affect the final sensitivity, proving this requires a technical study^{49} that is beyond the scope of this paper.

### APPENDIX J: GRAVITATIONAL WAVES ACTING ON A BOSE–EINSTEIN CONDENSATE

Consider a laboratory in a drag-free satellite, so that it is in free fall in the total gravitational field of the Earth and any GWs that may be present. In a sufficiently small region of space, we can choose coordinates $(t,x)$, such that the metric is flat, even in the presence of GWs. Expanding this metric to second order in $|xi|$, we have^{47}

where

and $g\mu \nu =\eta \mu \nu +h\mu \nu $, with $|h\mu \nu |\u226a1$. An Earthbound detector would not, in general, be in free fall with respect to the Earth's gravity, and so, in this case, the above metric must be adjusted as follows:^{47}

which is referred to as the proper detector frame,^{47,86} where *a* is the acceleration and $\Sigma $ is the rotation due to the Earth. However, a GW detector should be isolated as much as possible from the Earth's gravity and Newtonian noise. For example, the acceleration *a* can be compensated for by a suspension mechanism, and we can consider high enough GW frequencies that we can neglect Newtonian noise, so that only the GW contributes to the Riemann tensor at quadratic order. Then, we can use the freely falling frame metric (J1) to describe our experiment with only the GW entering the Riemann tensor.

In linearized theory, the Riemann tensor is invariant rather than just covariant, so we can compute it in any frame. It is simplest to choose the transverse-traceless (TT), where $h0\mu =0,\u2009hii=0,and\u2009\u2202jhij=0$, and taking a continuous monochromatic GW of frequency Ω propagating along the z-axis, *h _{ij}* is

Assuming that only the GW contributes through the Riemann tensor; we obtain the line element *ds _{D}* in our “lab” frame,

where $h+(t)=h+\u2009sin\u2009(\Omega t)\u2009and\u2009h\xd7(t)=h\xd7\u2009sin\u2009(\Omega t)$, and we have ignored the *z* dependence of $h\mu \nu TT$, given the large GW wavelengths compared to the size of the laboratory. Writing this metric as $g\mu \nu =:\eta \mu \nu \u2212hD\mu \nu $, we have

The Lagrangian density for a quantum, relativistic, massive, complex scalar field in a trap potential $v$, a spacetime with metric $g\mu \nu $, and $\Phi \u03024$ interaction potential is^{85}

Taking our above lab frame, with $g00=\u22121\u2212hD00,\u2009g0i=gi0=0,\u2009gij=\eta ij$, and $\u2212g=1\u221212hD00$, then

The BEC is non-relativistic limit, and so we should take $\Phi \u0302=\u210f2m\Psi \u0302e\u2212imc2t/\u210f$ and apply $|\u22020\Psi \u0302|\u226amc|\Psi \u0302|/\u210f?$ to obtain

where $g:=\u210f2\lambda /(2m)$. We must also take $c\u2192\u221e$, which, given that $h00\u221d1/c2$, reduces the above to

The Hamiltonian density is given by

where

The quantum field Hamiltonian for a non-relativistic Bose gas interacting with a GW can, therefore, be well-approximated by

where

which, for an effective one-dimensional BEC, reduces to

where we have used $h+(t)=h+\u2009sin\u2009(\Omega t)$ and defined $\epsilon :=h+$.

As noted in Appendix I, this result can also be derived by observing that a GW just acts as a Newtonian force obeying Hooke's law in the proper detector frame,^{47} with $Fi=mh\xa8ijTT\xi j/2$, where *ξ ^{j}* is the coordinate separation. The BEC will, thus, be subjected to such a force, and we can equivalently treat this as applying a quadratic potential to the BEC. The results using this potential coincide with the results in Ref. 18 obtained by solving the Klein–Gordon equation in the TT frame with an effective metric that depends on the GW spacetime metric, the BEC density, and speed of sound.

## References

**. See, e.g., Ref. 25, for an alternative convention.**

*x**q*and

*p*representations of the symplectic matrices corresponding to all the unitary processes involved in the interferometer.

^{56}The level of squeezing can then be calculated through $tanh(r)=e\u2212\beta /2$, where $\beta :=\u210f\omega max/(kBTBH)$, with

^{56}$TBH=0.36mcout2$ and $cs=0.57\u2009mm\u2009s\u22121$.

*et al.*, “