# Probing the quantum nature of gravity using a Bose-Einstein condensate
“Erste Abhandlung”

Soham Sen<EMAIL_ADDRESS>Department of Astrophysics and High Energy
Physics, S. N. Bose National Centre for Basic Sciences, JD Block, Sector-III,
Salt Lake City, Kolkata-700 106, India  Sunandan Gangopadhyay
<EMAIL_ADDRESS>Department of Astrophysics and High Energy
Physics, S. N. Bose National Centre for Basic Sciences, JD Block, Sector-III,
Salt Lake City, Kolkata-700 106, India

###### Abstract

The effect of noise induced by gravitons has been investigated using a Bose-
Einstein condensate. The general complex scalar field theory with a quadratic
self-interaction term has been considered in the presence of a gravitational
wave. The gravitational wave perturbation is then considerd as a sum of
discrete Fourier modes in the momentum space. Varying the action and making
use of the principle of least action, one obtains two equations of motion
corresponding to the gravitational perturbation and the time-dependent part of
the pseudo-Goldstone boson. Coming to an operatorial representation and
quantizing the phase space variables via appropriately introduced canonincal
commutation relations between the canonically conjugate variables
corresponding to the graviton and bosonic part of the total system, one
obtains a proper quantum gravity setup. Then we obtain the Bogoliubov
coefficients from the solution of the time-dependent part of the pseudo-
Goldstone boson and construct the covariance metric for the bosons initially
being in a squeezed state. The entries of the covariance matrix now involves a
stochastic contribution which results in an operatorial stochastic structure
of the quantum Fisher information. Using the stochastic average of the Fisher
information, we obtain a lower bound on the amplitude parameter of the
gravitational wave. As the entire calculation is done at zero temperature, the
bosonic system, by construction, will behave as a Bose-Einstein condensate.
For a Bose-Einstein condensate with a single mode, we observe that the lower
bound of the expectation value of the square of the uncertainty in the
amplitude measurement does not become infinite when the total observational
term approaches zero. It always has a finite value if the gravitons are
initially in a squeezed state with high enough squeezing. In order to sum over
all possible momentum modes, we next consider a noise term with a suitable
Gaussian weight factor which decays over time. We then obtain the lower bound
on the final expectation value of the square of the variance in the amplitude
parameter. Because of the noise induced by the graviton, there is a minimum
value of the measurement time below which it is impossible to detect any
gravitational wave using a Bose-Einstein condensate. Finally, we consider
interaction between the phonon modes of the Bose-Einstein condensate which
results in a decoherence. We observe that the decoherence effect becomes
significant for gravitons with minimal squeezing.

##  I Introduction

The derivation of the Planck’s radiation law by Satyendranath Bose in 1924 [1]
led to the introduction of the Bose statistics. Albert Einstein in this time
frame of 1924 and 1925 [2, 3, 4] extended this idea to matter systems which
led to the idea of a Bose gas governed by Bose statistics. Einstein also
proposed the existence of a new state of matter which was later termed as
Bose-Einstein condensate. The idea of a Bose-Einstein condensation is that if
a bosonic system (even a boson gas) is cooled below a critical temperature all
the bosons occupy the ground state energy level of the system. The matter
waves start superposing when the de-Broglie wavelength is larger than the
interatomic distance of the individual atoms and eventually at the moment of
crossing the critical temperature, all the matter waves superpose to form a
single wave function occupying the ground state of the bosonic system. This
phenomena is termed as Bose-Einstein condensation. Experimentally, Bose-
Einstein condensation was first detected in 1995 in a gas of Rubidium atoms
[5] and later in a gas of Sodium atoms in the same year [6]. Since then people
have tried to evolve the method for producing Bose-Einstein condensates and
use it for various physical applications. Another important aspect of
theoretical physics was the experimental detection of gravitational waves [7,
8, 9] which has opened up a new era of theoretical physics involving the
sculpturization of subatomic or lower physical length scales via the use of
gravitational waves. Gravitational wave detection by using atom interferometry
has been proposed quite a some time ago in [10, 11]. Recently in [12], a
gravitational wave detector using a Bose-Einstein condensate has been proposed
where a zero temperature quasi (1+1)-dimensional Bose-Einstein condensate with
fluctuating boundary conditions has been considered. An alternative
calculation considering the interaction of a nonrelativistic Bose-Einstein
condensate with a gravitational wave has been done in [13]. Later in [14, 15],
the idea propsed in [12] has been extended and made much more enhanced using a
$(3+1)$-dimensional zero temperature Bose-Einstein condensation and a decaying
gravitational wave template. The quantum Fisher information $H_{\varepsilon}$
was calculated by analyzing the fidelity between the individual squeezed
phonon states. The quantum Fisher information carries the amount of
information carried by the gravitational wave.

Recently, another aspect of high energy physics has emerged where the
stochastic effect of the noise of gravitons from a linearized quantum gravity
setup has been oberved [16, 17, 18, 19, 20, 21, 22, 23]. In [16, 17, 18], an
interferometer detector has been modelled by means of a freely falling pair of
particles in a slightly curved background. The perturbation over the Minkowski
spacetime has then been decomposed into its discrete Fourier modes in
$(1+1)$-dimension. Following a path integral approach the influence functional
of the gravitons over the detector system was calculated and varying the
action with respect to the detector degrees of freedom, the geodesic deviation
equation was obtained, which had the structure of a Langevin-like stochastic
differential equation. In these works, it has been shown that if a graviton is
initially in a squeezed state then it may be possible to detect signatures
from the detector-graviton interaction in future generation of gravitational
wave detectors. Another set of analysis were done in $(3+1)$-spacetime
dimensions and a canonical approach was followed. Similar stochastic Langevin-
like equations were obtained [19] and an indirect detection of gravitons by
means of decoherence was proposed in [20]. Similar but unique stochastic
effects has been observed in several other analyses [21, 22, 23].

The primary motivation of this work is to unveil the effects of the noise
induced by the gravitons on a homogeneous Bose-Einstein condensate in
$(3+1)$-spacetime dimensions at zero temperature. To carry out the analysis,
we need to start with the combined action comprising of the action describing
the Bose-Einstein condensate in curved spacetime and the Einstein-Hilbert
action. Here we have got rid of all the heavy fields in the theory as they
will have very small contribution towards the overall dynamics of the theory.
From [14, 15], we already know that a BEC is susceptible to gravitational wave
when the resonance condition is matched. As the gravitational fluctuation in
our analysis is now quantized, we expect to observe more subtle effects of the
gravitons on the phonons. Such small effects can lead to a BEC state which
will be incorporate such noise fluctuations. If one can now find a way to
trace such signatures of fluctuations due to graviton-BEC interaction, it will
suffice as an indirect detection of gravitons. In our analysis, we investigate
the BEC-graviton interaction using quantum metrological techniques and we
consider the quantum Fisher information to be the primary tool for indicating
quantum gravity signatures. Due to such noise fluctuations, one needs to look
at now the stochastic average of the quantum-gravity modified Fisher
information. The square root of the stochastic average of the Fisher
information will give the minimum value of the standard deviation in the
amplitude parameter of the gravitational wave. We investigate the form of the
quantum gravity modified Fisher information for squeezed graviton states which
also will highlight the primary analysis of the paper. Later we have
considered a scenario when the noise fluctuation is controlled by a Gaussian
decay factor. Next, we investigate on whether the BEC will be a good candidate
for extracting signatures of quantum gravity and compared it with required
sensitivity data from space based gravitational wave observatories. We have
finally investigated the effect of decoherence due to self interacting phonon
modes in the quantum gravity modified Fisher information.

Our paper is organized as follows. In section (II), we obtain the total action
of the system. In section (III), we discuss the noise induced by the gravitons
in the Bose-Einstein condensate. Later in section (IV) we discuss the quantum
metrology and obtain the quantum Fisher information for the system. We
consider a different noise template in section (V). In section (VI) we
consider decoherence due to interacting phonon modes and finally in section
VII, we summarize our results.

##  II Action of the system

In this section, we shall obtain the total action for the system in which a
gravitational wave is interacting with a self interacting scalar field theory
describing bosons. We work in the mostly positive signature for the metric.
The background metric can be thought of as a small perturbation on the flat
Minkowski background. The background metric is given by

$g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu}$ (1)

where $\eta_{\mu\nu}=diag\\{-1,1,1,1\\}$. If we consider the speed of light to
be unity, then the Einstein-Hilbert action can be written as

$S_{EH}=\frac{1}{16\pi G}\int d^{4}x\sqrt{-g}R$ (2)

with $R$ being the Ricci scalar and $g=\text{det}(g_{\mu\nu})$. Up to
quadratic order in the perturbation term in eq.(1), we can recast the Einstein
Hilbert action as follows

$\begin{split}S_{EH}\simeq&\frac{1}{64\pi G}\int d^{4}x~{}(h_{\mu\nu}\Box
h^{\mu\nu}-h\Box h+2h^{\mu\nu}\partial_{\mu}\partial_{\nu}h\\\
&-2h_{\mu\alpha}\partial_{\kappa}\partial^{\alpha}h^{\mu\kappa})~{}.\end{split}$
(3)

Now we shall make use of the gauge symmetry of the perturbation term given by

$h_{\mu\nu}=\bar{h}_{\mu\nu}+\partial_{\mu}\xi_{\nu}+\partial_{\nu}\xi_{\mu}~{}.$
(4)

Using this, we now impose the transverse-traceless gauge conditions given by

$\displaystyle\partial_{\kappa}\bar{h}^{\kappa\zeta}=0~{},\bar{h}^{\kappa}_{\kappa}=0~{},~{}k_{\rho}\bar{h}^{\rho\zeta}=0$
(5)

with $k_{\rho}=\delta_{\rho}^{0}$ being a constant time-like vector. In the
transverse-traceless gauge, the form of the Einstein-Hilbert action in eq.(3)
can be recast as

$S_{EH}=-\frac{1}{8\kappa^{2}}\int
d^{4}x~{}\partial_{\kappa}\bar{h}_{ij}\partial^{\kappa}\bar{h}^{ij}$ (6)

where $\kappa=\sqrt{8\pi G}$. One can now make use of a Fourier mode
decomposition of the fluctuation term $\bar{h}_{ij}$ inside a box of volume
$V$ as

$\bar{h}_{ij}(t,\mathbf{x})=\frac{2\kappa}{\sqrt{V}}\sum_{\mathbf{k},s}h^{s}(t,\mathbf{k})e^{i\mathbf{k}\cdot\mathbf{x}}\epsilon^{s}_{ij}(\mathbf{k})~{}.$
(7)

It is imperative to know that
$\bar{h}_{ij}(t,\mathbf{x})=\bar{h}_{ij}^{*}(t,\mathbf{x})$ as
$\bar{h}_{ij}(t,\mathbf{x})$ is a real quantity. Now making use of the Fourier
mode decomposition in eq.(7) and the reality condition of the fluctuation
term, we can recast the Einstein-Hilbert action in eq.(6) as

$\begin{split}S_{\text{EH}}=\frac{1}{2}\sum_{\textbf{k},s}\int
dt\left(\bigr{|}\dot{h}^{s}(t,\mathbf{k})\bigr{|}^{2}-k^{2}\bigr{|}h^{s}(t,\mathbf{k})\bigr{|}^{2}\right)~{}.\end{split}$
(8)

The Lagrangian density for a complex scalar bosonic field with a self
interaction term (in natural units) can be written as

$\begin{split}\mathcal{L}&=\nabla_{\mu}\phi^{\dagger}\nabla^{\mu}\phi+m^{2}\phi^{\dagger}\phi+\lambda\left(\phi^{\dagger}\phi\right)^{2}\\\
&=g^{\mu\nu}\partial_{\mu}\phi^{\dagger}(t,\mathbf{x})\partial_{\nu}\phi(t,\mathbf{x})+m^{2}|\phi(t,\mathbf{x})|^{2}+\lambda\left|\phi(t,\mathbf{x})\right|^{4}\end{split}$
(9)

where $m$ gives the mass of the bosonic field and
$\lambda\left|\phi\right|^{4}$ gives the self interaction term for the bosons.
Now eq.(9) effectively describes the Lagrangian density of a Bose-Einstein
condensate (BEC) as we are doing a zero-temperature field theory. Note that
the Lagrangian density presented in eq.(9) is a bit different from the one
presented in [14, 15] as we are working explicitly in a mostly positive
signature.

We now consider a homogeneous BEC and write $\phi$ as
$\phi(t,\mathbf{x})=e^{i\chi(t,\mathbf{x})}\varphi(t,\mathbf{x})$, where
$\chi(t,\mathbf{x})$ and $\varphi(t,\mathbf{x})$ are both real. Substituting
this relation back in eq.(9), we obtain the modified Lagrangian density as

$\begin{split}\mathcal{L}=g^{\mu\nu}\partial_{\mu}\varphi\partial_{\nu}\varphi+\varphi^{2}g^{\mu\nu}\partial_{\mu}\chi\partial_{\nu}\chi+m^{2}\varphi^{2}+\lambda\varphi^{4}~{}.\end{split}$
(10)

Here, $\varphi(t,\mathbf{x})$ is the heavy field, hence we extremize the
Lagrangian density in eq.(10) with respect to $\varphi$ as

$\begin{split}\frac{\partial\mathcal{L}}{\partial\varphi}=2\varphi
g^{\mu\nu}\partial_{\mu}\chi\partial_{\nu}\chi+2m^{2}\varphi+4\lambda\varphi^{3}&=0\\\
\implies
2\varphi\left(g^{\mu\nu}\partial_{\mu}\chi\partial_{\nu}\chi+m^{2}+2\lambda\varphi^{2}\right)&=0~{}.\end{split}$
(11)

As $\varphi$ is an arbitrary scalar field, it is possible to write down the
extremization condition from eq.(11) as

$\begin{split}\varphi^{2}=-\frac{1}{2\lambda}\left(g^{\mu\nu}\partial_{\mu}\chi\partial_{\nu}\chi+m^{2}\right)~{}.\end{split}$
(12)

Substituting the above relation back in the Lagrangian density in eq.(10) we
get,

$\begin{split}\mathcal{L}&=g^{\mu\nu}\partial_{\mu}\varphi\partial_{\nu}\varphi+\varphi^{2}\left(g^{\mu\nu}\partial_{\mu}\chi\partial_{\nu}\chi+m^{2}\right)+\lambda(\varphi^{2})^{2}\\\
&=g^{\mu\nu}\partial_{\mu}\varphi\partial_{\nu}\varphi-\frac{1}{4\lambda}\left(g^{\mu\nu}\partial_{\mu}\chi\partial_{\nu}\chi+m^{2}\right)^{2}~{}.\end{split}$
(13)

The primary focus of the analysis lies in the low frequency regime and as a
result the heavy field $\varphi$ can be integrated out from the theory [24,
14]. As a result, we can define a new Lagrangian density with an effective
minus sign as

$\mathcal{L}_{\text{BEC}}=\frac{1}{4\lambda}\left(g_{\mu\nu}\partial^{\mu}\chi\partial^{\nu}\chi+m^{2}\right)^{2}~{}~{}.$
(14)

Corresponding to this new Lagrangian, one can write down the total action of
the matter part of the system (which is the BEC coupled to the gravity) as

$S_{\text{BEC}}=\int d^{4}x\sqrt{-g}\mathcal{L}_{\text{BEC}}$ (15)

where $g=\text{det}[g_{\mu\nu}]$. If $\pi(t,\mathbf{x})\in\mathbb{R}$ denotes
the BEC phonons then in terms of these pseudo-Goldstone bosons, we can express
$\chi$ as

$\chi(t,\mathbf{x})=-\tilde{\sigma}t+\pi(t,\mathbf{x})=\tilde{\sigma}x_{\mu}\delta^{\mu}_{~{}0}+\pi(t,\mathbf{x})$
(16)

where $x^{0}=t$ and $x_{0}=g_{0\mu}x^{\mu}=-t$.

The background metric has the form given by

$g_{\mu\nu}=\begin{pmatrix}-1&0&0&0\\\ 0&1+h_{+}(t)&h_{\times}(t)&0\\\
0&h_{\times}(t)&1-h_{+}(t)&0\\\ 0&0&0&1\end{pmatrix}$ (17)

where $h_{+}$ and $h_{\times}$ denote the plus and cross polarizations of the
gravitational wave, propagating in the $z$ direction. From eq.(17), we obtain,
$\sqrt{-g}=\sqrt{1-(h_{+}^{2}+h_{\times}^{2})}\simeq
1+\mathcal{O}(h_{\mu\nu}h^{\mu\nu})$. Using the decomposition in eq.(16) and
using the expansion of $\sqrt{-g}$, we can recast the action in eq.(15) as

$\begin{split}S_{\text{BEC}}&\simeq\frac{1}{4\lambda}\int
d^{4}x\left(g_{\mu\nu}(\tilde{\sigma}\delta^{\mu}_{~{}0}+\partial^{\mu}\pi)(\tilde{\sigma}\delta^{\nu}_{~{}0}+\partial^{\nu}\pi)+m^{2}\right)^{2}\end{split}$
(18)

where we have kept $\sqrt{-g}$ upto the leading order term and neglected the
$\mathcal{O}(h_{\mu\nu}h^{\mu\nu})$ contribution. One can neglect the higher
order derivative terms [14] while expanding eq.(18). Hence, eq.(18) can be
recast in the following form

$\begin{split}S_{\text{BEC}}=&\int\frac{d^{4}x}{2\lambda}\bigr{[}(3\tilde{\sigma}^{2}-m^{2})\dot{\pi}^{2}-(\tilde{\sigma}^{2}-m^{2})g_{ij}\partial^{i}\pi\partial^{j}\pi\bigr{]}\\\
&+\int\frac{d^{4}x}{4\lambda}\left[4\tilde{\sigma}(\tilde{\sigma}^{2}-m^{2})\dot{\pi}+(\tilde{\sigma}^{2}-m^{2})^{2}\right]\\\
=&\int\frac{d^{4}x}{2\lambda}\bigr{[}(3\tilde{\sigma}^{2}-m^{2})\dot{\pi}^{2}-(\tilde{\sigma}^{2}-m^{2})g_{ij}\partial^{i}\pi\partial^{j}\pi\bigr{]}\\\
&+\int\frac{d^{3}x}{\lambda}\left[\tilde{\sigma}(\tilde{\sigma}^{2}-m^{2})\pi\right]\\\
\simeq&\int\frac{d^{4}x}{2\lambda}\bigr{[}(3\tilde{\sigma}^{2}-m^{2})\dot{\pi}^{2}-(\tilde{\sigma}^{2}-m^{2})g_{ij}\partial^{i}\pi\partial^{j}\pi\bigr{]}\end{split}$
(19)

where in the second line of the above equation we have got rid of the non-
dynamical contributions and in the final line we have made use of the fact
that $\pi(t,\mathbf{x})$ vanishes at the boundary. We now make use of the
ansatz for the pseudo-Goldstone boson

$\pi(t,\mathbf{x})=\sum_{\mathbf{k}_{\beta}}e^{i\mathbf{k}_{\beta}\cdot\mathbf{x}}\psi_{\mathbf{k}_{\beta}}(t)~{}.$
(20)

As $\pi(t,\mathbf{x})\in\mathbb{R}$, we know that
$\pi(t,\mathbf{x})=\pi^{*}(t,\mathbf{x})$. This reality condition leads us to
the relation
$\sum_{\mathbf{k}_{\beta}}e^{i\mathbf{k}_{\beta}\cdot\mathbf{x}}\psi_{\mathbf{k}_{\beta}}(t)=\sum_{\mathbf{k}_{\beta}}e^{i\mathbf{k}_{\beta}\cdot\mathbf{x}}\psi^{*}_{-\mathbf{k}_{\beta}}(t).$
The above relation implies
$\psi_{\mathbf{k}_{\beta}}(t)=\psi^{*}_{-\mathbf{k}_{\beta}}(t)$
$\forall~{}\mathbf{k}_{\beta}$. Throughout the analysis, we have neglected the
spatial dependence of the gravitational fluctuation (this assumption has also
been adopted in the classical treatment of [14]). Using the above condition
and the Fourier mode decomposition of the gravitational fluctuation term from
eq.(7), eq.(19) can be recast as

$\begin{split}S_{\text{BEC}}\simeq&\frac{1}{2\lambda}\int
d^{4}x\left[(3\tilde{\sigma}^{2}-m^{2})\dot{\pi}(t,\mathbf{x})\dot{\pi}^{*}(t,\mathbf{x})-(\tilde{\sigma}^{2}-m^{2})\left(\eta_{ij}+\bar{h}_{ij}(t,0)\right)\partial^{i}\pi(t,\mathbf{x})\partial^{j}\pi^{*}(t,\mathbf{x})\right]\\\
=&\frac{1}{2\lambda}\int
d^{4}x\biggr{[}(3\tilde{\sigma}^{2}-m^{2})\sum\limits_{\mathbf{k}_{\beta},\mathbf{k}^{\prime}_{\beta}}e^{i(\mathbf{k}_{\beta}-\mathbf{k}^{\prime}_{\beta})\cdot\mathbf{x}}\dot{\psi}_{\mathbf{k}_{\beta}}(t)\dot{\psi}^{*}_{\mathbf{k}^{\prime}_{\beta}}(t)-(\tilde{\sigma}^{2}-m^{2})\Bigr{[}\eta_{ij}+\frac{2\kappa}{\sqrt{V}}\sum\limits_{\mathbf{k},s}h_{\mathbf{k},s}(t)\epsilon^{s}_{ij}(\mathbf{k})\Bigr{]}\\\
&\times\sum\limits_{\mathbf{k}_{\beta},\mathbf{k}^{\prime}_{\beta}}(ik_{\beta}^{i})(-i{k^{\prime}}_{\beta}^{j})e^{i(\mathbf{k}_{\beta}-\mathbf{k}^{\prime}_{\beta})\cdot\mathbf{x}}\psi_{\mathbf{k}_{\beta}}(t)\psi^{*}_{\mathbf{k}^{\prime}_{\beta}}(t)\biggr{]}\\\
=&\frac{1}{2\lambda}\int
dt\biggr{[}(3\tilde{\sigma}^{2}-m^{2})\sum\limits_{\mathbf{k}_{\beta},\mathbf{k}^{\prime}_{\beta}}\dot{\psi}_{\mathbf{k}_{\beta}}(t)\dot{\psi}^{*}_{\mathbf{k}^{\prime}_{\beta}}(t)\int
d^{3}x~{}e^{i(\mathbf{k}_{\beta}-\mathbf{k}^{\prime}_{\beta})\cdot\mathbf{x}}-(\tilde{\sigma}^{2}-m^{2})\Bigr{[}\eta_{ij}+\frac{2\kappa}{\sqrt{V}}\sum\limits_{\mathbf{k},s}h_{\mathbf{k},s}(t)\epsilon^{s}_{ij}(\mathbf{k})\Bigr{]}\\\
&\times\sum\limits_{\mathbf{k}_{\beta},\mathbf{k}^{\prime}_{\beta}}(ik_{\beta}^{i})(-i{k^{\prime}}_{\beta}^{j})\psi_{\mathbf{k}_{\beta}}(t)\psi^{*}_{\mathbf{k}^{\prime}_{\beta}}(t)\int
d^{3}x~{}e^{i(\mathbf{k}_{\beta}-\mathbf{k}^{\prime}_{\beta})\cdot\mathbf{x}}\biggr{]}\\\
=&\frac{1}{2\lambda}\int
dt\biggr{[}(3\tilde{\sigma}^{2}-m^{2})\sum\limits_{\mathbf{k}_{\beta},\mathbf{k}^{\prime}_{\beta}}\dot{\psi}_{\mathbf{k}_{\beta}}(t)\dot{\psi}^{*}_{\mathbf{k}^{\prime}_{\beta}}(t)V_{\beta}\delta_{\mathbf{k}_{\beta},\mathbf{k}^{\prime}_{\beta}}-(\tilde{\sigma}^{2}-m^{2})\Bigr{[}\eta_{ij}+\frac{2\kappa}{\sqrt{V}}\sum\limits_{\mathbf{k},s}h_{\mathbf{k},s}(t)\epsilon^{s}_{ij}(\mathbf{k})\Bigr{]}\\\
&\times\sum\limits_{\mathbf{k}_{\beta},\mathbf{k}^{\prime}_{\beta}}k_{\beta}^{i}{k^{\prime}}_{\beta}^{j}\psi_{\mathbf{k}_{\beta}}(t)\psi^{*}_{\mathbf{k}^{\prime}_{\beta}}(t)V_{\beta}\delta_{\mathbf{k}_{\beta},\mathbf{k}^{\prime}_{\beta}}\biggr{]}\\\
\implies S_{\text{BEC}}=&\frac{V_{\beta}}{2\lambda}\int
dt\biggr{[}(3\tilde{\sigma}^{2}-m^{2})\sum\limits_{\mathbf{k}_{\beta}}\bigr{|}\dot{\psi}_{\mathbf{k}_{\beta}}(t)\bigr{|}^{2}-(\tilde{\sigma}^{2}-m^{2})\Bigr{[}\eta_{ij}+\frac{2\kappa}{\sqrt{V}}\sum\limits_{\mathbf{k},s}h_{\mathbf{k},s}(t)\epsilon^{s}_{ij}(\mathbf{k})\Bigr{]}\sum\limits_{\mathbf{k}_{\beta}}k_{\beta}^{i}{k}_{\beta}^{j}\bigr{|}\psi_{\mathbf{k}_{\beta}}(t)\bigr{|}^{2}\biggr{]}\end{split}$
(21)

where in the penultimate line we have made use of the normalization condition
for the pseudo-Goldstone bosons inside a box of volume $V_{\beta}$ as $\int
d^{3}x~{}e^{i(\mathbf{k}_{\beta}-\mathbf{k}^{\prime}_{\beta})\cdot\mathbf{x}}=V_{\beta}\delta_{\mathbf{k}_{\beta},\mathbf{k}^{\prime}_{\beta}}$.
Here, $\delta_{\mathbf{k}_{\beta},\mathbf{k}^{\prime}_{\beta}}$ is the
abbreviated form of
$\delta_{\mathbf{k}_{\beta},\mathbf{k}^{\prime}_{\beta}}\equiv\delta_{{k}^{1}_{\beta},{k^{\prime}}^{1}_{\beta}}\times\delta_{{k}^{2}_{\beta},{k^{\prime}}^{2}_{\beta}}\times\delta_{{k}^{3}_{\beta},{k^{\prime}}^{3}_{\beta}}$.
Such a box normalization of the pseudo-Goldstone bosons is quite intuitive in
a sense that experimentally a Bose-Einstein condensate is formed in a very
confined region (e.g. making use of harmonic trap potentials) which can
resemble the shape of a small box. From the dispersion relation of the BEC
[14], we know that $\omega_{\beta}^{2}\simeq c_{s}^{2}k_{\beta}^{2}$, where
$c_{s}^{2}=\frac{\tilde{\sigma}^{2}-m^{2}}{3\tilde{\sigma}^{2}-m^{2}}$
(denoting the square of the speed of sound) and $k_{\beta}\ll m$ with
$c_{s}\ll 1$ (in natural units, $c=1$). We can further recast eq.(21) as

$\begin{split}S_{\text{BEC}}=&\gamma_{\beta}\int
dt\biggr{[}\sum\limits_{\mathbf{k}_{\beta}}\bigr{|}\dot{\psi}_{\mathbf{k}_{\beta}}(t)\bigr{|}^{2}-c_{s}^{2}\Bigr{[}\eta_{ij}+\frac{2\kappa}{\sqrt{V}}\sum\limits_{\mathbf{k},s}h_{\mathbf{k},s}(t)\epsilon^{s}_{ij}(\mathbf{k})\Bigr{]}\sum\limits_{\mathbf{k}_{\beta}}k_{\beta}^{i}k_{\beta}^{j}\bigr{|}\psi_{\mathbf{k}_{\beta}}(t)\bigr{|}^{2}\biggr{]}\end{split}$
(22)

where $\gamma_{\beta}$ has the dimension of length in natural units and
$\gamma_{\beta}\equiv\frac{V_{\beta}}{2\lambda}(3\tilde{\sigma}^{2}-m^{2})$.
Combining the Einstein-Hilbert action from eq.(8) with that of the action for
the BEC from eq.(22), we obtain the total action of the system as

$\begin{split}S=&S_{\text{EH}}+S_{\text{BEC}}\\\
=&\frac{1}{2}\sum_{\textbf{k},s}\int
dt\left(\bigr{|}\dot{h}^{s}(t,\mathbf{k})\bigr{|}^{2}-k^{2}\bigr{|}h^{s}(t,\mathbf{k})\bigr{|}^{2}\right)+\gamma_{\beta}\int
dt\sum\limits_{\mathbf{k}_{\beta}}\biggr{[}\bigr{|}\dot{\psi}_{\mathbf{k}_{\beta}}(t)\bigr{|}^{2}-c_{s}^{2}\Bigr{[}\eta_{ij}+\frac{2\kappa}{\sqrt{V}}\sum\limits_{\mathbf{k},s}h_{\mathbf{k},s}(t)\epsilon^{s}_{ij}(\mathbf{k})\Bigr{]}k_{\beta}^{i}k_{\beta}^{j}\bigr{|}\psi_{\mathbf{k}_{\beta}}(t)\bigr{|}^{2}\biggr{]}~{}.\end{split}$
(23)

We start by varying the action given in eq.(23) in terms of the complex
conjugate of the time dependent part of the pseudo-Goldstone boson
corresponding to individual momentum modes and the complex conjugate of the
individual Fourier mode of the graviton. Using the principle of least action
($\frac{\delta S}{\delta\psi^{*}_{\mathbf{k}_{\beta}}}=0$), we obtain a
dynamical equation or the equation of motion corresponding to the time
dependent part of the pseudo-Goldstone boson as

$\begin{split}\ddot{\psi}_{\mathbf{k}_{\beta}}(t)+c_{s}^{2}\Bigr{[}\eta_{ij}+\frac{2\kappa}{\sqrt{V}}\sum\limits_{\mathbf{k},s}h_{\mathbf{k},s}(t)\epsilon^{s}_{ij}(\mathbf{k})\Bigr{]}{k}^{i}_{\beta}{k}^{j}_{\beta}\psi_{\mathbf{k}_{\beta}}(t)&=0~{}.\end{split}$
(24)

Similarly, by using the principle of least action after extremixing the action
with respect to the variable $h^{*}_{\mathbf{k},s}$, we get

$\begin{split}\ddot{h}_{\mathbf{k},s}(t)+k^{2}h_{\mathbf{k},s}(t)=-\frac{4\gamma_{\beta}\kappa
c_{s}^{2}}{\sqrt{V}}\epsilon^{s*}_{ij}(\mathbf{k})\sum\limits_{\mathbf{k}_{\beta}}k_{\beta}^{i}k_{\beta}^{j}\left\lvert\psi_{\mathbf{k}_{\beta}}(t)\right\rvert^{2}\end{split}$
(25)

where we have made use of the reality condition of $\bar{h}_{ij}(t,0)$. With
the two equations of motion in hand, we can now move towards writing down a
quantum mechanical model of the gravitational wave-BEC system.

##  III Graviton induced noise in the BEC

In this section, our primary aim is to quantize the theory. The simplest way
of quantizing the gravitational part is to raise the Fourier modes of the
spacetime fluctuation $\bar{h}_{ij}$ to operator status and impose a suitable
canonical commutation relation among $\hat{h}_{\mathbf{k},s}(t)$ and its
canonically conjugate variable in the phase space. For the bosonic part, it
can be a bit tricky. One can just raise the time dependent part of the pseudo-
Goldstone bosons to operator status and impose suitable commutation relation
among the two canonically conjugate variables in the phase space. The second
way is to quantize it in the momentum space and raise the momentum variables
to operator status with the suitable use of canonical commutation relations.
In this first paper of our set of two works, we shall make use of the first
procedure, where we only quantize $\psi_{\mathbf{k}_{\beta}}(t)$ in its
entirety and impose commutation relation between
$\hat{\psi}_{\mathbf{k}_{\beta}}(t)$ and its canonically conjugate variable.
The independent Fourier mode of the spacetime fluctuation can be decomposed
into two parts. One is the classical contribution which is obtained by taking
the expectation value of the mode operator and a quantum fluctuation term. It
is now possible to write down the fluctuation term of the mode operator
corresponding to a momentum value $\mathbf{k}$ in the interaction picture as
[19]

$\delta\hat{h}^{I}_{\mathbf{k},s}(t)=\hat{h}_{\mathbf{k},s}^{I}(t)-h_{\text{cl}}^{s}(\mathbf{k},t)$
(26)

with the classical component given as
$h_{\text{cl}}^{s}(\mathbf{k},t)=\bigr{\langle}\hat{h}_{\mathbf{k},s}^{I}(t)\bigr{\rangle}$
where the expectation is taken with respect to the initial state of the
graviton. In principle the classical contribution has non-vanishing
contribution if the graviton is initially in coherent, squeezed vacuum,
thermal, and similar other combinatorial states. Here,
$\delta\hat{h}^{I}_{\mathbf{k},s}(t)$ is the gravitational quantum fluctuation
[19]. The quantum field in the interaction picture in terms of the creation
and annihilation operators can be represented as

$\hat{h}_{\mathbf{k},s}^{I}(t)=u_{k}(t)\hat{a}_{s}(\mathbf{k})+u_{k}^{*}(t)\hat{a}^{\dagger}_{s}(-\mathbf{k})$
(27)

where $k=|\mathbf{k}|$ with the mode function, $u_{k}(t)$, satisfying the
following normalization condition

$-iu_{k}(t)\overset{\leftrightarrow}{\partial}_{t}u_{k}^{*}(t)=-i\left(u_{k}(t)\dot{u}_{k}^{*}(t)-\dot{u}_{k}(t)u_{k}^{*}(t)\right)=1~{}.$
(28)

In case of Minkowski vacuum, the creation and annihilation operators in
eq.(27), satisfy the following commutation relation

$\begin{split}[\hat{a}_{s}(\mathbf{k}),\hat{a}^{\dagger}_{s^{\prime}}(\mathbf{k}^{\prime})]&=\delta_{s,s^{\prime}}\delta_{\mathbf{k},\mathbf{k}^{\prime}},\\\
[\hat{a}_{s}(\mathbf{k}),\hat{a}_{s^{\prime}}(\mathbf{k}^{\prime})]&=[\hat{a}^{\dagger}_{s}(\mathbf{k}),\hat{a}^{\dagger}_{s^{\prime}}(\mathbf{k}^{\prime})]=0~{}.\end{split}$
(29)

Raising $h_{\mathbf{k},s}(t)$ and $\hat{\psi}_{\mathbf{k}_{\beta}}(t)$ to
operator status, we can recast eq.(25) in a quantum mechanical representation
as

$\begin{split}\ddot{\hat{h}}_{\mathbf{k},s}(t)+k^{2}\hat{h}_{\mathbf{k},s}(t)=-\frac{4\gamma_{\beta}\kappa
c_{s}^{2}}{\sqrt{V}}\epsilon^{s*}_{ij}(\mathbf{k})\sum\limits_{\mathbf{k}_{\beta}}k_{\beta}^{i}k_{\beta}^{j}\left\lvert\hat{\psi}_{\mathbf{k}_{\beta}}(t)\right\rvert^{2}\end{split}$
(30)

where
$\left\lvert\hat{\psi}_{\mathbf{k}_{\beta}}(t)\right\rvert^{2}=\hat{\psi}_{\mathbf{k}_{\beta}}^{\dagger}(t)\hat{\psi}_{\mathbf{k}_{\beta}}(t)$.
Making use of the Green’s function technique, it is possible to write down the
solution of the above equation as

$\begin{split}\hat{h}_{\mathbf{k},s}(t)=&\hat{h}^{I}_{\mathbf{k},s}(t)-\frac{4\gamma_{\beta}\kappa
c_{s}^{2}}{\sqrt{V}}\epsilon^{s*}_{ij}(\mathbf{k})\int_{0}^{t}dt^{\prime}\frac{\sin(k(t-t^{\prime}))}{k}\sum\limits_{\mathbf{k}^{\prime}_{\beta}}{{k}^{\prime}_{\beta}}^{i}{{k}^{\prime}_{\beta}}^{j}\bigr{|}\hat{\psi}_{\mathbf{k}^{\prime}_{\beta}}(t^{\prime})\bigr{|}^{2}\\\
=&h_{\text{cl}}^{s}(\mathbf{k},t)+\delta\hat{h}^{I}_{\mathbf{k},s}(t)-\frac{4\gamma_{\beta}\kappa
c_{s}^{2}}{\sqrt{V}}\epsilon^{s*}_{ij}(\mathbf{k})\sum\limits_{\mathbf{k}^{\prime}_{\beta}}{k^{\prime}_{\beta}}^{i}{k^{\prime}_{\beta}}^{j}\int_{0}^{t}dt^{\prime}\frac{\sin(k(t-t^{\prime}))}{k}\bigr{|}\hat{\psi}_{\mathbf{k}^{\prime}_{\beta}}(t^{\prime})\bigr{|}^{2}~{}.\end{split}$
(31)

We can also write down the quantum mechanical version of eq.(24) as

$\begin{split}\ddot{\hat{\psi}}_{\mathbf{k}_{\beta}}(t)+c_{s}^{2}\Bigr{[}\eta_{ij}+\frac{2\kappa}{\sqrt{V}}\sum\limits_{\mathbf{k},s}\hat{h}_{\mathbf{k},s}(t)\epsilon^{s}_{ij}(\mathbf{k})\Bigr{]}{k}^{i}_{\beta}{k}^{j}_{\beta}\hat{\psi}_{\mathbf{k}_{\beta}}(t)&=0\end{split}$
(32)

where in the last line of the above equation, we have made use of eq.(26). One
can now regulate the mode summations, corresponding to the gravitational wave
part, via the use of an ultraviolet (UV) cut-off. Substituting the form of
$\hat{h}_{\mathbf{k},s}(t)$ from eq.(31), we can write down eq.(32) as

$\begin{split}&\ddot{\hat{\psi}}_{\mathbf{k}_{\beta}}(t)+c_{s}^{2}\eta_{ij}k_{\beta}^{i}k_{\beta}^{j}\hat{\psi}_{\mathbf{k}_{\beta}}(t)+c_{s}^{2}\biggr{(}\frac{2\kappa}{\sqrt{V}}\sum\limits_{s}\smashoperator[]{\sum_{\begin{subarray}{c}{\mathbf{k}}\\\
{|\mathbf{k}|\leq\Omega_{m}}\end{subarray}}^{}}h_{\text{cl}}^{s}(\mathbf{k},t)\epsilon^{s}_{ij}(\mathbf{k})\biggr{)}k_{\beta}^{i}k_{\beta}^{j}\hat{\psi}_{\mathbf{k}_{\beta}}(t)+c_{s}^{2}\biggr{(}\frac{2\kappa}{\sqrt{V}}\sum\limits_{s}\smashoperator[]{\sum_{\begin{subarray}{c}{\mathbf{k}}\\\
{|\mathbf{k}|\leq\Omega_{m}}\end{subarray}}^{}}\delta\hat{h}^{I}_{\mathbf{k},s}(t)\epsilon^{s}_{ij}(\mathbf{k})\biggr{)}k_{\beta}^{i}k_{\beta}^{j}\hat{\psi}_{\mathbf{k}_{\beta}}(t)\\\
&-\frac{8\gamma_{\beta}\kappa^{2}c_{s}^{4}}{V}\sum\limits_{s}\smashoperator[]{\sum_{\begin{subarray}{c}{\mathbf{k}}\\\
{|\mathbf{k}|\leq\Omega_{m}}\end{subarray}}^{}}\epsilon^{s*}_{ij}(\mathbf{k})\epsilon^{s}_{lm}(\mathbf{k})\biggr{(}\int_{0}^{t}dt^{\prime}\frac{\sin(k(t-t^{\prime}))}{k}\sum\limits_{\mathbf{k}^{\prime}_{\beta}}{k^{\prime}_{\beta}}^{i}{k^{\prime}_{\beta}}^{j}\bigr{|}\hat{\psi}_{\mathbf{k}^{\prime}_{\beta}}(t^{\prime})\bigr{|}^{2}\biggr{)}k_{\beta}^{l}k_{\beta}^{m}\hat{\psi}_{\mathbf{k}_{\beta}}(t)=0~{}.\end{split}$
(33)

We will henceforth use the UV-regulated mode summations for the gravitational
wave part throughout this work. From an experimental scenario this is quite
logical as an usual gravitational wave detector can detect frequencies upto a
maximum value.

One can now define two new quantities as

$\displaystyle h_{ij}^{\text{cl}}(t,\mathbf{x})\equiv$
$\displaystyle\frac{2\kappa}{\sqrt{V}}\sum\limits_{s}\smashoperator[]{\sum_{\begin{subarray}{c}{\mathbf{k}}\\\
{|\mathbf{k}|\leq\Omega_{m}}\end{subarray}}^{}}h_{\text{cl}}^{s}(\mathbf{k},t)e^{i\mathbf{k}\cdot\mathbf{x}}\epsilon^{s}_{ij}(\mathbf{k})$
(34) $\displaystyle\delta\hat{N}_{ij}(t)\equiv$
$\displaystyle\frac{2\kappa}{\sqrt{V}}\sum\limits_{s}\smashoperator[]{\sum_{\begin{subarray}{c}{\mathbf{k}}\\\
{|\mathbf{k}|\leq\Omega_{m}}\end{subarray}}^{}}\delta\hat{h}^{I}_{\mathbf{k},s}(t)\epsilon^{s}_{ij}(\mathbf{k})~{}.$
(35)

One can now introduce the projection tensors to proceed further in this
analysis as follows

$\begin{split}\mathcal{P}_{ij}=\delta_{ij}-\frac{k_{i}k_{j}}{k^{2}}\end{split}$
(36)

and making use of them, one can write down the following relation involving
the sum over all polarizations of the product of two polarization tensors

$\begin{split}\sum\limits_{s}\epsilon^{s*}_{ij}(\mathbf{k})\epsilon^{s}_{lm}(\mathbf{k})=\frac{1}{2}\left[\mathcal{P}_{il}\mathcal{P}_{jm}+\mathcal{P}_{im}\mathcal{P}_{jl}-\mathcal{P}_{ij}\mathcal{P}_{lm}\right]~{}.\end{split}$
(37)

Using eq.(s)(34,35) and eq.(37), we can recast eq.(33) as

$\begin{split}&\ddot{\hat{\psi}}_{\mathbf{k}_{\beta}}(t)+c_{s}^{2}\left(\eta_{ij}+h_{ij}^{\text{cl}}(t,0)+\delta\hat{N}_{ij}(t)\right)k_{\beta}^{i}k_{\beta}^{j}\hat{\psi}_{\mathbf{k}_{\beta}}(t)\\\
&-\frac{\xi_{\beta}}{V}\smashoperator[]{\sum_{\begin{subarray}{c}{\mathbf{k}}\\\
{|\mathbf{k}|\leq\Omega_{m}}\end{subarray}}^{}}\left(\mathcal{P}_{il}\mathcal{P}_{jm}+\mathcal{P}_{im}\mathcal{P}_{jl}-\mathcal{P}_{ij}\mathcal{P}_{lm}\right)\sum\limits_{\mathbf{k}^{\prime}_{\beta}}{k^{\prime}_{\beta}}^{i}{k^{\prime}_{\beta}}^{j}\\\
&\times\left(\int_{0}^{t}dt^{\prime}\frac{\sin(k(t-t^{\prime}))}{k}\bigr{|}\hat{\psi}_{\mathbf{k}^{\prime}_{\beta}}(t^{\prime})\bigr{|}^{2}\right)k_{\beta}^{l}k_{\beta}^{m}\hat{\psi}_{\mathbf{k}_{\beta}}(t)=0\end{split}$
(38)

where $\xi_{\beta}=4\gamma_{\beta}\kappa^{2}c_{s}^{4}$, and
$\frac{\xi_{\beta}}{V}$ is a dimensionless number. The summation over the
graviton modes can be converted in a continuous mode integral as
$\frac{1}{V}\smashoperator[]{\sum_{{\mathbf{k}};{|\mathbf{k}|\leq\Omega_{m}}}^{}}\rightarrow\frac{1}{(2\pi)^{3}}\int^{\Omega_{m}}d^{3}k$.
In this integral $d^{3}k$ can be converted to corresponding spherical
coordinates in Fourier space as $k^{2}dk\sin\theta d\theta
d\phi=k^{2}dkd\Omega$. The angular integrals are given by [19]

$\begin{split}&\int d\Omega=4\pi,~{}\int
d\Omega~{}k^{i}k^{j}=\frac{4\pi}{3}\delta_{ij},\\\ &\int
d\Omega~{}k^{i}k^{j}k^{l}k^{m}=\frac{4\pi}{15}\left(\delta^{ij}\delta^{lm}+\delta^{il}\delta^{jm}+\delta^{im}\delta^{jl}\right).\end{split}$
(39)

Using the angular integrals and making use of the discrete to continuous mode
conversion rule, we obtain the following result for the summation terms
involving the projection tensors,
$\sum_{\mathbf{k}}\left(\mathcal{P}_{il}\mathcal{P}_{jm}+\mathcal{P}_{im}\mathcal{P}_{jl}-\mathcal{P}_{ij}\mathcal{P}_{lm}\right)$,
to be

$\begin{split}&\frac{V}{(2\pi)^{3}}\int^{\Omega_{m}}_{0}dkk^{2}\int
d\Omega\left[\mathcal{P}_{il}\mathcal{P}_{jm}+\mathcal{P}_{im}\mathcal{P}_{jl}-\mathcal{P}_{ij}\mathcal{P}_{lm}\right]=\frac{V}{5\pi^{2}}\left(\int_{0}^{\Omega_{m}}dkk^{2}\right)\left(\delta_{il}\delta_{jm}+\delta_{im}\delta_{jl}-\frac{2}{3}\delta_{ij}\delta_{lm}\right)~{}.\end{split}$
(40)

Making use of eq.(40) and performing the $k$ integral, one can recast eq.(38)
as

$\begin{split}&\ddot{\hat{\psi}}_{\mathbf{k}_{\beta}}(t)+c_{s}^{2}\left(\eta_{ij}+h_{ij}^{\text{cl}}(t,0)+\delta\hat{N}_{ij}(t)\right)k_{\beta}^{i}k_{\beta}^{j}\hat{\psi}_{\mathbf{k}_{\beta}}(t)-\frac{\xi_{\beta}}{5\pi^{2}}\left(\delta_{il}\delta_{jm}+\delta_{im}\delta_{jl}-\frac{2}{3}\delta_{ij}\delta_{lm}\right)\sum\limits_{\mathbf{k}^{\prime}_{\beta}}{k^{\prime}_{\beta}}^{i}{k^{\prime}_{\beta}}^{j}\\\
&\times\int_{0}^{t}dt^{\prime}\left(\frac{\sin(\Omega_{m}(t-t^{\prime}))}{(t-t^{\prime})^{2}}-\frac{\Omega_{m}\cos(\Omega_{m}(t-t^{\prime}))}{(t-t^{\prime})}\right)\bigr{|}\hat{\psi}_{\mathbf{k}^{\prime}_{\beta}}(t^{\prime})\bigr{|}^{2}k_{\beta}^{l}k_{\beta}^{m}\hat{\psi}_{\mathbf{k}_{\beta}}(t)=0~{}.\end{split}$
(41)

Absorbing the Kronecker-deltas, we can resimplify the above equation as

$\begin{split}&\ddot{\hat{\psi}}_{\mathbf{k}_{\beta}}(t)+c_{s}^{2}\left(\eta_{ij}+h_{ij}^{\text{cl}}(t,0)+\delta\hat{N}_{ij}(t)\right)k_{\beta}^{i}k_{\beta}^{j}\hat{\psi}_{\mathbf{k}_{\beta}}(t)-\frac{2\xi_{\beta}}{5\pi^{2}}\int_{0}^{t}dt^{\prime}\left(\frac{\sin(\Omega_{m}(t-t^{\prime}))}{(t-t^{\prime})^{2}}-\frac{\Omega_{m}\cos(\Omega_{m}(t-t^{\prime}))}{(t-t^{\prime})}\right)\\\
&\times\sum_{\mathbf{k}^{\prime}_{\beta}}\left((\mathbf{k}_{\beta}\cdot\mathbf{k}^{\prime}_{\beta})^{2}-\frac{1}{3}k_{\beta}^{2}{k^{\prime}}^{2}_{\beta}\right)\bigr{|}\hat{\psi}_{\mathbf{k}^{\prime}_{\beta}}(t^{\prime})\bigr{|}^{2}\hat{\psi}_{\mathbf{k}_{\beta}}(t)=0~{}.\end{split}$
(42)

Our next aim is to obtain a solution for $\hat{\psi}_{\mathbf{k}_{\beta}}(t)$
via solving the above equation. The time-dependent part of the pseudo-
Goldstone boson can be divided into three parts,
$\hat{\psi}_{\mathbf{k}_{\beta}}(t)=\psi^{(0)}_{\mathbf{k}_{\beta}}(t)+\psi^{h_{\text{cl}}}_{\mathbf{k}_{\beta}}(t)+\hat{\psi}^{(1)}_{\mathbf{k}_{\beta}}(t)$.
Here, $\psi_{\mathbf{k}_{\beta}}^{(0)}(t)$ denotes the unperturbed classical
part of the solution and $\psi_{\mathbf{k}_{\beta}}^{h_{\text{cl}}}(t)$
denotes the first order solution corresponding to the classical gravitational
perturbation. The final part $\hat{\psi}^{(1)}_{\mathbf{k}_{\beta}}(t)$
encodes the solution corresponding to the quantum fluctuations of the
gravitons. It is important to note that the decomposition of the solution is
done in a way such that the operatorial contribution can be separated from the
classical part.

The zeroth order classical equation from eq.(42) reads

$\begin{split}\ddot{\psi}^{(0)}_{\mathbf{k}_{\beta}}(t)+c_{s}^{2}k_{\beta}^{2}\psi^{(0)}_{\mathbf{k}_{\beta}}(t)=0\end{split}$
(43)

which has a solution of the form

$\begin{split}\psi^{(0)}_{\mathbf{k}_{\beta}}(t)&=\mathcal{a}e^{-ic_{s}k_{\beta}t}+\mathcal{b}e^{ic_{s}k_{\beta}t}=\mathcal{a}e^{-i\omega_{\beta}t}+\mathcal{b}e^{i\omega_{\beta}t}~{}.\end{split}$
(44)

It is now quite intuitive to get rid of the negative energy modes and set
$\mathcal{b}=0$, and as a result of the normalization condition, we get
$\mathcal{a}=1$. The first order classical equation of motion from eq.(42) can
be written as

$\begin{split}\ddot{\psi}^{h_{\text{cl}}}_{\mathbf{k}_{\beta}}(t)+c_{s}^{2}k_{\beta}^{2}{\psi}^{h_{\text{cl}}}_{\mathbf{k}_{\beta}}(t)=-c_{s}^{2}h^{\text{cl}}_{ij}(t,0)k_{\beta}^{i}k_{\beta}^{j}\psi^{(0)}_{\mathbf{k}_{\beta}}(t)+\mathcal{f}_{\beta}(t)\end{split}$
(45)

where $\mathcal{f}_{\beta}(t)$ is given by

$\begin{split}\mathcal{f}_{\beta}(t)&=\frac{2\xi_{\beta}}{5\pi^{2}}\sum_{\mathbf{k}^{\prime}_{\beta}}\left[(\mathbf{k}_{\beta}\cdot\mathbf{k}^{\prime}_{\beta})^{2}-\frac{1}{3}k_{\beta}^{2}{k^{\prime}}^{2}_{\beta}\right]\psi^{(0)}_{\mathbf{k}_{\beta}}(t)\int_{0}^{t}dt^{\prime}\\\
\times&\biggr{[}\frac{\sin(\Omega_{m}(t-t^{\prime}))}{(t-t^{\prime})^{2}}-\frac{\Omega_{m}\cos(\Omega_{m}(t-t^{\prime}))}{(t-t^{\prime})}\biggr{]}\bigr{|}{\psi}^{(0)}_{\mathbf{k}^{\prime}_{\beta}}(t^{\prime})\bigr{|}^{2}\\\
&=\frac{2\xi_{\beta}}{5\pi^{2}}\left[\Omega_{m}-\frac{\sin(\Omega_{m}t)}{t}\right]\sum_{\mathbf{k}^{\prime}_{\beta}}\Bigr{[}(\mathbf{k}_{\beta}\cdot\mathbf{k}^{\prime}_{\beta})^{2}\\\
&-\frac{1}{3}k_{\beta}^{2}{k^{\prime}}^{2}_{\beta}\Bigr{]}e^{-i\omega_{\beta}t}~{}.\end{split}$
(46)

To proceed further and in order to make the analysis simpler, we restrict
ourselves to plus polarization of the gravitational wave only. As a result, we
already know that $\epsilon^{\times}_{ij}(\mathbf{k})=0$ ($\forall
i,j=\\{1,2,3\\}$),
$\epsilon^{+}_{11}(\mathbf{k})=-\epsilon^{+}_{22}(\mathbf{k})$,
$\epsilon^{+}_{33}(\mathbf{k})=0$, and $\epsilon^{+}_{ij}(\mathbf{k})=0$
$\forall i\neq j$. Making use of eq.(34), we can therefore write
$h^{\text{cl}}_{11}(t,0)=-h^{\text{cl}}_{22}(t,0)=h^{\text{cl}}(t,0)$. Hence,
eq.(45) can be recast as

$\begin{split}\ddot{\psi}^{h_{\text{cl}}}_{\mathbf{k}_{\beta}}(t)+\omega_{\beta}^{2}\psi^{h_{\text{cl}}}_{\mathbf{k}_{\beta}}(t)&=-\omega_{\beta}^{2}\mathcal{k}_{0}^{2}h^{\text{cl}}(t,0)e^{-i\omega_{\beta}t}+\mathcal{f}_{\beta}(t)\end{split}$
(47)

where
$\mathcal{k}_{0}^{2}=\frac{{k_{\beta}}_{x}^{2}-{k_{\beta}}_{y}^{2}}{k_{\beta}^{2}}$.
Corresponding to the Green’s function equation
$(\frac{d^{2}}{dt^{2}}+\omega_{\beta}^{2})\mathcal{G}(t-t^{\prime})=\delta(t-t)^{\prime}$,
we obtain the analytical form of the Green’s function to be

$\mathcal{G}(t-t^{\prime})=\frac{1}{\omega_{\beta}}\sin(\omega_{\beta}(t-t^{\prime}))\Theta(t-t^{\prime})$
(48)

with $\Theta(t-t^{\prime})$ denoting the Heaviside theta function. The
solution of eq.(45) can then be obtained as

$\begin{split}&\psi^{h_{\text{cl}}}_{\mathbf{k}_{\beta}}(t)=\mathcal{a}_{h}e^{-i\omega_{\beta}t}+\mathcal{b}_{h}e^{i\omega_{\beta}t}-\mathcal{k}_{0}^{2}\int_{-\infty}^{t}dt^{\prime}\omega_{\beta}h^{\text{cl}}(t^{\prime},0)\\\
&\times
e^{-i\omega_{\beta}t^{\prime}}\sin(\omega_{\beta}[t-t^{\prime}])+\int_{-\infty}^{t}dt^{\prime}\frac{\sin(\omega_{\beta}[t-t^{\prime}])}{\omega_{\beta}}\mathcal{f}_{\beta}(t^{\prime})~{}.\end{split}$
(49)

The simplest choice for the undetermined constants is to set
$\mathcal{a}_{h}=\mathcal{b}_{h}=0$. As before, we can also get
$\delta\hat{N}_{11}(t)=-\delta\hat{N}_{22}(t)=\delta\hat{N}(t)$ and other
components of the $\delta\hat{N}_{ij}(t)$ tensor is zero. The final dynamical
equation, involving the operators only (the quantum-gravitational time
evolution equation), reads

$\begin{split}\ddot{\hat{\psi}}^{(1)}_{\mathbf{k}_{\beta}}(t)+\omega_{\beta}^{2}\hat{\psi}^{(1)}_{\mathbf{k}_{\beta}}(t)\simeq-\omega_{\beta}^{2}\mathcal{k}_{0}^{2}\delta\hat{N}(t).\end{split}$
(50)

As $\hat{\psi}^{(1)}_{\mathbf{k}_{\beta}}(t)$ is an operator corresponding to
quantum gravity consideration, the operatorial contribution from
$\mathcal{f}_{\beta}(t)$ becomes way smaller compared to the other terms in
eq.(50). Again setting the random constants to zero and using the Green’s
function technique, we arrive at the solution of eq.(50) as follows

$\begin{split}\hat{\psi}^{(1)}_{\mathbf{k}_{\beta}}(t)=-\mathcal{k}_{0}^{2}\int_{-\infty}^{t}dt^{\prime}\omega_{\beta}\sin(\omega_{\beta}(t-t^{\prime}))e^{-i\omega_{\beta}t^{\prime}}\delta\hat{N}(t^{\prime})~{}.\end{split}$
(51)

As the quantum fluctuations purely arises because of the interaction of the
BEC with the gravitons, it is safe to assume $\delta\hat{N}(t)=0$, $\forall
t<0$. As a result, eq.(51) reduces to an integral whose limits are from $0$ to
$t$. Combining eq.(s)(44,47) and eq.(51) and using the specific values for the
constant, we obtain the complete solution for the time dependent part of the
pseudo-Goldstone boson as

$\hat{\psi}_{\mathbf{k}_{\beta}}(t)=e^{-i\omega_{\beta}t}-\omega_{\beta}\mathcal{k}_{0}^{2}\int^{t}_{-\infty}dt^{\prime}e^{-i\omega_{\beta}t^{\prime}}\sin(\omega_{\beta}(t-t^{\prime}))\left(h^{\text{cl}}(t^{\prime},0)+\delta\hat{N}(t^{\prime})\right)+\frac{1}{\omega_{\beta}}\int^{t}_{-\infty}dt^{\prime}\sin(\omega_{\beta}(t-t^{\prime}))\mathcal{f}_{\beta}(t^{\prime})~{}.$
(52)

This is one of the pivotal results in our paper and it signifies the fact that
the time dependent part of the pseudo-Goldstone boson now explicitly depends
upon a quantum fluctuation parameter which is a direct consequence of the
interaction of the supercooled BEC with gravitons. The subsequent discussion
in this paper will be based upon this important result that we have obtained.
We shall now proceed to calculate the form of the above solution when specific
incoming gravitational wave template is used and obtain the corresponding
Bogoliubov coefficients. At first, it is important to note that in the last
term in eq.(52) the integral can be assumed to operate in a very small time
limit as the interaction of BEC with a gravitational wave will occur for a
very small amount of time. From eq.(46), we observe that if $t$ is very small
then
$\Omega_{m}-\frac{\sin(\Omega_{m}t)}{t}\simeq\Omega_{m}-\frac{\Omega_{m}t}{t}=0$.
Hence, for a simpler analysis, we indeed can set $\mathcal{f}_{\beta}(t)$
equals to zero in eq.(52). The classical template of the gravitational wave
can be used as $h^{\text{cl}}(t^{\prime},0)=\varepsilon
e^{-\frac{t^{2}}{\tau^{2}}}\cos(\Omega t)$, where $\tau$ indicates the
duration of time for capturing a single measurement of the gravitational wave.
From eq.(52), we replace $h^{\text{cl}}(t^{\prime},0)$ by the above analytical
form and obtain
$I_{h}(t)\equiv\varepsilon\int_{-\infty}^{t}dt^{\prime}e^{-i\omega_{\beta}t^{\prime}}\sin(\omega_{\beta}(t-t^{\prime}))e^{-\frac{{t^{\prime}}^{2}}{\tau^{2}}}\cos(\Omega
t^{\prime})$. We shall make the upper limit to $\infty$ and argue that it is a
quite good approximation for a such a gravitational wave template and we plot
the two behaviour where the upper limit of integration is finite in one case
and infinite in the other case in Figure (1). The quantum gravitational
fluctuation term in eq.(52) is very small in amplitude. Hence, we can simply
assume that the noise overall can be modelled by that of the final value of
the fluctuation at time $t$. Hence, we use the relation that $\int
dt^{\prime}\delta\hat{N}(t^{\prime})g(t^{\prime})\simeq\int
dt^{\prime}\delta\hat{N}(t)g(t^{\prime})$ with $g(t)$ being a function of
time. In a simple experimental setup, if the measurement time corresponding to
a single BEC mode is $t=\mathcal{t}$ then the solution in eq.(52) can be
recast in the final form as

Figure 1: $I_{h}(t)$ is plotted for $\omega_{\beta}=1\text{ Hz}$,
$\Omega=2\text{ Hz}$ and $\tau=10$ sec when the upper limit of integration is
finite and infinite as well.

$\begin{split}\hat{\psi}_{\mathbf{k}_{\beta}}(t)&=\left(1-\frac{\mathcal{k}_{0}^{2}}{4}(1+2i\omega_{\beta}\mathcal{t})\delta\hat{N}(\mathcal{t})\right)e^{-i\omega_{\beta}t}+\left(\frac{\varepsilon\mathcal{k}_{0}^{2}}{4}\sqrt{\pi}\omega_{\beta}\tau\left(e^{-\frac{\tau^{2}}{4}(\Omega-2\omega_{\beta})^{2}}-e^{-\frac{\tau^{2}}{4}(\Omega+2\omega_{\beta})^{2}}\right)+\frac{\mathcal{k_{0}^{2}}}{4}\delta\hat{N}(\mathcal{t})\right)e^{i\omega_{\beta}t}\\\
&=\hat{\alpha}^{\beta}(\mathcal{t})e^{-i\omega_{\beta}t}+\hat{\beta}^{\beta}(\mathcal{t})e^{i\omega_{\beta}t}\end{split}$
(53)

where the coefficients $\hat{\alpha}^{\beta}(\mathcal{t})$ and
$\hat{\beta}^{\beta}(\mathcal{t})$ are defined as

$\displaystyle\hat{\alpha}^{\beta}(\mathcal{t})$
$\displaystyle\equiv\alpha^{\beta}+\delta\hat{\alpha}^{\beta}(\mathcal{t})=1-\frac{{k_{\beta}}_{x}^{2}-{k_{\beta}}_{y}^{2}}{4k_{\beta}^{2}}(1+2i\omega_{\beta}\mathcal{t})\delta\hat{N}(\mathcal{t})$
(54) $\displaystyle\hat{\beta}^{\beta}(\mathcal{t})$
$\displaystyle\equiv\beta^{\beta}+\delta\hat{\beta}^{\beta}(\mathcal{t})=\frac{{k_{\beta}}_{x}^{2}-{k_{\beta}}_{y}^{2}}{4k_{\beta}^{2}}\sqrt{\pi}\varepsilon\omega_{\beta}\tau\left(e^{-\frac{\tau^{2}}{4}(\Omega-2\omega_{\beta})^{2}}-e^{-\frac{\tau^{2}}{4}(\Omega+2\omega_{\beta})^{2}}\right)+\frac{{k_{\beta}}_{x}^{2}-{k_{\beta}}_{y}^{2}}{4k_{\beta}^{2}}\delta\hat{N}(\mathcal{t})~{}.$
(55)

It is quite natural to set $\mathcal{t}=\tau$ later on but for the moment we
leave them as two distinct numbers. It is straight forward to write down
$\delta\hat{\alpha}^{\beta}=\mathcal{c}^{\alpha}(\mathcal{t})\delta\hat{N}(\mathcal{t})$
and
$\delta\hat{\beta}^{\beta}=\mathcal{c}^{\beta}(\mathcal{t})\delta\hat{N}(\mathcal{t})$
such that
$\mathcal{c}^{\alpha}(\mathcal{t})=-\frac{\tilde{\varepsilon}}{4\epsilon}(1+2i\omega_{\beta}t)$,
and $\mathcal{c}^{\beta}(\mathcal{t})=\frac{\varepsilon}{4\varepsilon}$. Here,
a new quantity is defined,
$\tilde{\varepsilon}\equiv\varepsilon\frac{{k_{\beta}}_{x}^{2}-{k_{\beta}}_{y}^{2}}{k_{\beta}^{2}}$.
As the background spacetime is curved, it is evident from eq.(53) that
$\hat{\alpha}_{\beta}(\mathcal{t})$ and $\hat{\beta}_{\beta}(\mathcal{t})$ are
the Bogoliubov coefficients.

Eq.(s)(54,55) signifies that the Bogoliubov coefficients can now be decomposed
into two parts, one is the classical part which is same as that obtained in
[14], and the another part is a fluctuation term which is indeed purely a
quantum gravitational term. It is quite natural that as the background
spacetime fluctuations are quantized, the Bogoliubov coefficients will not be
constant numbers anymore rather the noise fluctuations will be embedded into
them. This is a very important observation in this paper. Our aim is to obtain
the variance in the measurement of the parameter $\varepsilon$. For the next
part of our analysis, we shall make use of quantum metrology techniques to
extract signatures of the noise induced by the gravitons on the BEC. It is
important to note that the metrology techniques followed in [14], will be a
bit different in this scenario as the Bogoliubov coefficients are now having a
fluctuation term. In order to inspect the quantum Fisher information term in
this scenario [14], we finally will need to take the stochastic average of the
same.

##  IV Quantum metrology and the noise of gravitons

###  IV.1 Calculating the covariance matrix for the single mode of an n mode
bosonic system

The method that we shall follow in this literature is the covariance matrix
formalism. In this subsection, we have given a pedagogical derivation of the
covariance matrix of a squeezed one mode Bose-Einstein condensate at zero
temperature. For a system of $n$ bosons, the position and the conjugate
momenta in terms of the ladder operators are given as

$\hat{x}^{\beta}_{j}=\sqrt{\frac{\hbar}{2m\omega_{\beta}}}(\hat{a}^{\beta}_{j}+\hat{a}^{\beta\dagger}_{j})~{},~{}~{}\hat{p}^{\beta}_{k}=i\sqrt{\frac{m\hbar\omega_{\beta}}{2}}(\hat{a}^{\beta\dagger}_{k}-\hat{a}^{\beta}_{k})~{}.$
(56)

If one imposes the commutation relation
$[\hat{a}^{\beta}_{k},\hat{a}^{\beta\dagger}_{k^{\prime}}]=\delta_{k,k^{\prime}}$,
it is quite straight forward to check that
$[\hat{x}^{\beta}_{j},\hat{p}^{\beta}_{k}]=i\hbar\delta_{jk}$ with
$j,k=1,\cdots,n$. Here, $\hat{a}_{k}^{\beta}$ denotes the annihilation
operator corresponding to the $k$-th energy state such that
$\hat{a}_{k}^{\beta}|m_{k}^{\beta}\rangle=\sqrt{m_{k}^{\beta}}|m_{k}^{\beta}-1\rangle$,
and $\hat{a}_{k}^{\beta\dagger}$ denotes the creation operator such that
$\hat{a}_{k}^{\beta\dagger}|m_{k}^{\beta}\rangle=\sqrt{m_{k}^{\beta}+1}|m_{k}^{\beta}+1\rangle$.
The vacuum state is defined as $\hat{a}_{k}^{\beta}|0^{\beta}_{k}\rangle=0$,
$\forall k$. One can now define a column vector of the form
$\mathcal{R}=\left(\mathcal{u}\hat{x}^{\beta}_{1},\frac{1}{\mathcal{v}}\hat{p}^{\beta}_{1},\cdots,\mathcal{u}\hat{x}^{\beta}_{n},\frac{1}{\mathcal{v}}\hat{p}^{\beta}_{n}\right)^{T}$
where $\mathcal{u}\equiv\sqrt{\frac{m\omega_{\beta}}{\hbar}}$ and
$\mathcal{v}=\sqrt{m\hbar\omega_{\beta}}$. It is quite straight forward to
check that
$[\mathcal{R},\mathcal{R}^{T}]\equiv\mathcal{R}\mathcal{R}^{T}-(\mathcal{R}\mathcal{R}^{T})^{T}=i\mathcal{O}$,
where $\mathcal{O}=\bigoplus\limits_{j=1}^{n}i\sigma_{2}$ with $\sigma_{2}$
denoting the second Pauli spin matrix and $\bigoplus$ denoting the direct sum
corresponding to the $n$-modes of the $n$-mode bosonic system. Here,
$\mathcal{O}$ is $2n\times 2n$ dimensional matrix. The covariant matrix
$\Sigma$ in terms of the $\mathcal{R}$ column matrix reads

$\Sigma_{ij}=\frac{1}{2}\langle\\{\mathcal{R}_{i},\mathcal{R}_{j}\\}\rangle-\langle\mathcal{R}_{i}\rangle\langle\mathcal{R}_{i}\rangle$
(57)

where the expectation value is taken with resect to the density matrix of the
$n$ mode bosons. For a single mode, we can write down
$\mathcal{R}^{k}=\left(\mathcal{u}\hat{x}^{\beta}_{k},\frac{1}{\mathcal{v}}\hat{p}_{k}^{\beta}\right)^{T}$
where only the $k$-th mode is being considered. It is straight forward to
check that for a single-mode vacuum state of bosons in thermal equilibrium,
the density matrix reads ($k$-th bosonic mode is considered)

$\begin{split}\hat{\rho}^{k}&=\frac{e^{-\text{\ss}\hat{H}}}{\text{tr}[e^{-\text{\ss}\hat{H}}]}\\\
&=\frac{\sum\limits_{n_{k}^{\beta}=0}^{\infty}|n_{k}^{\beta}\rangle\langle
n_{k}^{\beta}|e^{-\text{\ss}\hat{a}_{k}^{\beta\dagger}\hat{a}_{k}^{\beta}}}{\sum\limits_{m_{k}^{\beta}=0}^{\infty}\langle
m_{k}^{\beta}|e^{-\text{\ss}\hat{a}_{k}^{\beta\dagger}\hat{a}_{k}^{\beta}}|m_{k}^{\beta}\rangle}\\\
&=(1-e^{-\text{\ss}})\sum\limits_{n_{k}^{\beta}=0}^{\infty}|n_{k}^{\beta}\rangle\langle
n_{k}^{\beta}|e^{-\text{\ss}n_{k}^{\beta}}\\\
&=\frac{1}{1+\mathcal{N}}\sum\limits_{n_{k}^{\beta}=0}^{\infty}\left(\frac{\mathcal{N}}{1+\mathcal{N}}\right)^{m_{k}^{\beta}}|m_{k}^{\beta}\rangle\langle
m_{k}^{\beta}|\end{split}$ (58)

where $\text{\ss}=\frac{1}{k_{B}T}$ and
$\mathcal{N}=\frac{1}{e^{\text{\ss}}-1}$ with $k_{B}$ denoting the Boltzmann
constant, and $T$ denoting the equilibrium temperature. Using eq.(58), it is
straight forward to show that $\langle\mathcal{R}^{k}\rangle=0$. It is
therefore sufficient to calculate
$\frac{1}{2}\text{tr}\left[\\{\mathcal{R}^{k},{\mathcal{R}^{k}}^{T}\\}\hat{\rho}^{k}\right]$
in order to obtain the covariance matrix $\Sigma^{k}$ corresponding to the
$k$-th bosonic mode of the $n$-mode bosonic system. The anti-commutator
$\\{\mathcal{R}^{k},{\mathcal{R}^{k}}^{T}\\}$ reads

$\\{\mathcal{R}^{k},{\mathcal{R}^{k}}^{T}\\}=\begin{bmatrix}\frac{2m\omega_{\beta}}{\hbar}(\hat{x}^{\beta}_{k})^{2}&&\frac{1}{\hbar}\left(\hat{x}^{\beta}_{k}\hat{p}^{\beta}_{k}+\hat{p}^{\beta}_{k}\hat{x}^{\beta}_{k}\right)\\\
\frac{1}{\hbar}\left(\hat{p}^{\beta}_{k}\hat{x}^{\beta}_{k}+\hat{x}^{\beta}_{k}\hat{p}^{\beta}_{k}\right)&&\frac{2}{m\hbar\omega_{\beta}}(\hat{p}_{k}^{\beta})^{2}\end{bmatrix}.$
(59)

Using the above matrix structure, one can obtain the final form of the
covariance matrix corresponding to a single mode of an $n$-mode bosonic system
as

$\begin{split}\Sigma^{k}[\mathcal{N}(T)]&=\frac{1}{2}\text{tr}\left[\\{\mathcal{R}^{k},{\mathcal{R}^{k}}^{T}\\}\hat{\rho}^{k}\right]=\frac{2\mathcal{N}+1}{2}\begin{bmatrix}1&&0\\\
0&&1\end{bmatrix}~{}.\end{split}$ (60)

In the zero temperature limit, $\mathcal{N}=0$. If the entire system is at
zero temperature then effectively it denotes a BEC corresponding to the $k$-th
mode of the $n$-mode bosonic system. The covariance matrix corresponding to a
BEC from eq.(60) can be obtained using the zero temperature limit as

$\begin{split}\Sigma^{k}[\mathcal{N}(0)]&=\Sigma^{k}[0]\\\
&=\lim\limits_{T\rightarrow 0}\Sigma^{k}[\mathcal{N}(T)]\\\
&=\begin{bmatrix}\frac{1}{2}&&0\\\ 0&&\frac{1}{2}\end{bmatrix}~{}.\end{split}$
(61)

In order to enhance the feedback of the BEC from the gravitational wave, the
general idea is to squeeze the single mode bosons in the zero temperature
limit. For the next part of our analysis, we shall drop the $k$ superscript.
Under a squeezing by a parameter $r_{\text{sq.}}=re^{i\varphi}$, one can
obtain the following two relations

$\begin{split}\hat{S}(r)\hat{a}^{\beta}\hat{S}^{\dagger}(r)&=\hat{a}^{\beta}\cosh
r+\hat{a}^{\beta\dagger}e^{i\varphi}\sinh r~{},\\\
\hat{S}(r)\hat{a}^{\beta\dagger}\hat{S}^{\dagger}(r)&=\hat{a}^{\beta\dagger}\cosh
r+\hat{a}^{\beta}e^{-i\varphi}\sinh r~{}.\end{split}$ (62)

Applying the transformations to the vector $\mathcal{R}$, we obtain

$\begin{split}&\hat{S}(r)\mathcal{R}\hat{S}^{\dagger}(r)=\begin{pmatrix}\mathcal{u}\hat{S}(r)\hat{x}^{\beta}\hat{S}^{\dagger}(r)\\\
\frac{1}{\mathcal{v}}\hat{S}(r)\hat{p}^{\beta}\hat{S}^{\dagger}(r)\end{pmatrix}\\\
&=\begin{bmatrix}\cosh r+\cos\varphi\sinh r&&\sin\varphi\sinh r\\\
\sin\varphi\sinh r&&\cosh r-\cos\varphi\sinh
r\end{bmatrix}\begin{pmatrix}\mathcal{u}\hat{x}^{\beta}\\\
\frac{1}{\mathcal{v}}\hat{p}^{\beta}\end{pmatrix}\\\
&=\Xi_{\text{sq.}}(r)\mathcal{R}\end{split}$ (63)

where $\Xi_{\text{sq.}}(r)$ denotes the squeezing matrix. Hence, the
covariance matrix corresponding to the single-mode squeezed phonons of the BEC
reads 111Similar result for the squeezed covariance matrix has been reproduced
earlier in [14, 25] but a slightly different result was produced where the
off-diagonal elements of the matrix comes with a negative sign.

$\begin{split}&\Sigma_{\text{sq.}}[0]\\\
&=\Xi_{\text{sq.}}(r)~{}\Sigma[0]~{}\Xi^{T}_{\text{sq.}}(r)\\\
&=\frac{1}{2}\begin{bmatrix}\cosh 2r+\cos\varphi\sinh 2r&&\sin\varphi\sinh
2r\\\ \sin\varphi\sinh 2r&&\cosh 2r-\cos\varphi\sinh
2r\end{bmatrix}~{}.\end{split}$ (64)

When a gravitational wave interacts with the squeezed BEC, it will transform
the covariance matrix obtained in eq.(64) as [26]

$\begin{split}\tilde{\Sigma}_{k}(\tilde{\varepsilon})=\mathcal{M}_{kk}(\tilde{\varepsilon})\Sigma_{\text{sq.}}[0]\mathcal{M}_{kk}^{T}(\tilde{\varepsilon})+\sum\limits_{j\neq
k}\mathcal{M}_{kj}(\tilde{\varepsilon})\mathcal{M}^{T}_{kj}(\tilde{\varepsilon})\end{split}$
(65)

where $\tilde{\varepsilon}=\varepsilon\mathcal{k}_{0}^{2}$ and
$\mathcal{M}_{kj}(\tilde{\varepsilon})$ is given as [14, 26]

$\mathcal{M}_{kj}(\tilde{\varepsilon})=\begin{bmatrix}\Re[{\alpha}^{\beta}_{kj}-{\beta}^{\beta}_{kj}]&&\Im[{\alpha}^{\beta}_{kj}+{\beta}^{\beta}_{kj}]\\\
-\Im[{\alpha}^{\beta}_{kj}-{\beta}^{\beta}_{kj}]&&\Re[{\alpha}^{\beta}_{kj}+{\beta}^{\beta}_{kj}]\end{bmatrix}$
(66)

with ${\alpha}^{\beta}_{kj}$ and $\beta^{\beta}_{kj}$ denoting the classical
Bogoliubov coefficients. In our current analysis, the Bogoliubov coefficients
are operators involving a small fluctuation term. The expectation value of the
fluctuation term vanishes and the two point correlator has a non-vanishing
contribution. Hence, the individual elements of the matrix
$\mathcal{M}_{kj}(\varepsilon)$ will have an additional contribution from the
noise fluctuations. The modified symplectic matrix including the effects from
the noise fluctuation takes the form

$\begin{split}&\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\tilde{\mathcal{M}}}}}_{kj}(\tilde{\varepsilon})\equiv\mathcal{M}_{kj}(\tilde{\varepsilon})+\delta\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\mathcal{M}}}}_{kj}(\tilde{\varepsilon})\\\
&=\mathcal{M}_{kj}(\tilde{\varepsilon})+\begin{bmatrix}\Re[\mathcal{c}^{\alpha}_{kj}-\mathcal{c}^{\beta}_{kj}]&&\Im[\mathcal{c}^{\alpha}_{kj}+\mathcal{c}^{\beta}_{kj}]\\\
-\Im[\mathcal{c}^{\alpha}_{kj}-\mathcal{c}^{\beta}_{kj}]&&\Re[\mathcal{c}^{\alpha}_{kj}+\mathcal{c}^{\beta}_{kj}]\end{bmatrix}\delta\hat{N}(\mathcal{t})\end{split}$
(67)

where we have defined a new symbol
$\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{A}}}$ which indicates a
matrix $A$ with operators as its elements. It is essential to note that
$\delta\hat{N}(\mathcal{t})$ is a stochastic parameter, as a result we cannot
define its eigenvalues and the non-vanishing contribution comes only from the
two-point correlator of the stochastic term with the contribution being always
a real number. As a result, we can call it a stochastic operator and
$\delta\hat{M}_{kj}(\tilde{\varepsilon})$ carries the entire essence of the
stochastic operator. It is though interesting to note that
$\delta\hat{N}(\mathcal{t})$ has no well defined adjoint operator, as a result
it is complementary to consider it as real operator. We can now rewrite
eq.(65) in terms of the modified symplectic matrices with elements including
operators as

$\begin{split}\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\tilde{\Sigma}}}}_{k}(\tilde{\varepsilon})=\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\tilde{\mathcal{M}}}}}_{kk}(\tilde{\varepsilon})\Sigma_{\text{sq.}}[0]\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\tilde{\mathcal{M}}}}}_{kk}^{T}(\tilde{\varepsilon})+\sum\limits_{j\neq
k}\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\tilde{\mathcal{M}}}}}_{kj}(\tilde{\varepsilon})\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\tilde{\mathcal{M}}}}}^{T}_{kj}(\tilde{\varepsilon})~{}.\end{split}$
(68)

The Bogoliubov coefficients do not involve two different modes corresponding
to the $n$-mode bosonic system. Hence, it is straight forward to express the
two coefficients as
$\hat{\alpha}_{kj}^{\beta}=\delta_{kj}\hat{\alpha}^{\beta}$ and
$\hat{\beta}_{kj}^{\beta}=\delta_{kj}\hat{\beta}^{\beta}$. With the analytical
form of
$\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\tilde{\mathcal{M}}}}}_{kj}(\tilde{\varepsilon})$,
we are now in a position to calculate the error in measurement of the
parameter $\tilde{\varepsilon}$ in the next subsection.

###  IV.2 Quantum Fisher information

In a general scenario, classical measurement may suffice in precise
determination of parameters. In cases where quantum mechanical effects are
mostly in action, it may not be possible to precisely determine the outcome of
small parameter without using quantum mechanical measurement techniques. This
is known as quantum metrology.

####  IV.2.1 Cramér-Rao bound and the quantum Fisher information

Here, we give a brief introduction to the Cramér-Rao bound involving the
classical Fisher information and the techniques used to obtain the quantum
Fisher information. Consider a generalized measurement by a set of Hermitian
operators $\hat{\mathcal{G}}(\zeta)$ which are nonnegative and $\int
d\zeta~{}\hat{\mathcal{G}}(\zeta)=\hat{\mathbb{1}}$. If the probability
density for obtaining the result $\zeta$, when a parameter $\vartheta$ is
given, is
$p(\zeta|\vartheta)=\text{tr}[\hat{\mathcal{G}}(\zeta)\hat{\rho}(\vartheta)]$
then the classical Fisher information is defined by

$\begin{split}\mathcal{I}_{\vartheta}&\equiv\int
d\zeta~{}p(\zeta|\vartheta)\left[\frac{\partial\ln
p(\zeta|\vartheta)}{\partial\vartheta}\right]^{2}\\\
&=\int\frac{d\zeta}{p(\zeta|\vartheta)}\left[\frac{\partial
p(\zeta|\vartheta)}{\partial\vartheta}\right]^{2}~{}.\end{split}$ (69)

The minimum value in the error of the estimation of the parameter $\vartheta$
from the $\mathfrak{N}$ number of independent measurements, with the set of
results $\\{\zeta_{1},\zeta_{2},\cdots,\zeta_{\mathfrak{N}}\\}$, is obtained
using the Cramér-Rao bound to be [27]

$\langle(\Delta\vartheta)^{2}\rangle\geq\frac{1}{\mathfrak{N}\mathcal{I}_{\vartheta}}~{}.$
(70)

If one now considers $\vartheta$ to be a parameter corresponding to a quantum-
mechanical system, then the generalized form of the classical Fisher
information reads [27]

$\mathcal{I}_{\vartheta}=\int
d\zeta\frac{1}{\text{tr}[\hat{\mathcal{G}}(\zeta)\hat{\rho}(\vartheta)]}\text{tr}\left[\hat{\mathcal{G}}(\zeta)\frac{\partial\hat{\rho}(\vartheta)}{\partial\vartheta}\right]^{2}~{}.$
(71)

The quantum Fisher information, (when considering all measurements
$\\{\hat{\mathcal{G}}(\zeta)\\}$) reads [27]

$\mathcal{H}_{\vartheta}=\max\limits_{\\{\hat{\mathcal{G}}(\zeta)\\}}\mathcal{I}_{\vartheta}~{}.$
(72)

Hence, the maximum amount of information, one can extract after $\mathfrak{N}$
measurements is determined by the quantum Fisher information as

$\langle(\Delta\vartheta)^{2}\rangle\geq\frac{1}{\mathfrak{N}\mathcal{I}_{\vartheta}}\geq\frac{1}{\mathfrak{N}\mathcal{H}_{\vartheta}}~{}.$
(73)

For two states $\rho_{1}$ and $\rho_{2}$, the overlap between them is
determined by the fidelity
$\mathcal{F}(\rho_{1},\rho_{2})=\left(\text{tr}\left[\sqrt{\sqrt{\rho_{1}}\rho_{2}\sqrt{\rho}_{2}}\right]\right)$.
One can express the quantum Fisher information in eq.(72), in terms of the
Fidelity between two nearby states $\rho_{\vartheta}$ and
$\rho_{\vartheta+d\vartheta}$ as [26]

$\mathcal{H}_{\vartheta}=\frac{8(1-\sqrt{\mathcal{F}(\rho_{\vartheta},\rho_{\vartheta+d\vartheta})})}{d\vartheta^{2}}~{}.$
(74)

For Gaussian states, it is easier to use the covariance matrix approach than
the density matrix approach. Now, the overlap between two covariance matrices
$\Sigma_{1}$ and $\Sigma_{2}$ for a single mode bosonic systems reads [28]

$\mathcal{F}(\Sigma_{1},\Sigma_{2})=\frac{1}{\sqrt{\Lambda+\Delta}-\sqrt{\Lambda}}$
(75)

where

$\displaystyle\Lambda$
$\displaystyle=\frac{1}{4}\det\left[\Sigma_{1}+\frac{i}{2}\mathcal{O}\right]\det\left[\Sigma_{2}+\frac{i}{2}\mathcal{O}\right]~{},$
(76) $\displaystyle\Delta$
$\displaystyle=\frac{1}{4}\det\left[\Sigma_{1}+\Sigma_{2}\right]~{}.$ (77)

If $\vartheta$ is a very small parameter then it is possible to perturbatively
expand $\Sigma$ along with $\mathcal{F}$ and $\mathcal{H}_{\vartheta}$. We
briefly discuss the methodology presented in [26]. For a perturbative
calculation, the initial assumption is that the Bogoliubov coefficients can be
expanded upto second order in $\vartheta$ as

$\begin{split}\alpha_{ij}(\vartheta)&\simeq\alpha_{ij}^{(0)}+\vartheta\alpha_{ij}^{(1)}+\vartheta^{2}\alpha_{ij}^{(2)}\\\
\beta_{ij}(\vartheta)&\simeq\vartheta\beta_{ij}^{(1)}+\vartheta^{2}\beta_{ij}^{(2)}~{}.\end{split}$
(78)

The above expansion is applicable for any symplectic matrix $\mathcal{M}$ that
operates on $\Sigma$ to change it to some different matrix, such that the
uncertainty relation still holds true. As both the first and second order
moments of $\mathcal{R}$ can be expanded in this manner, one can express the
covariance matrix $\Sigma(\vartheta)$ upto order $\vartheta^{2}$ as

$\Sigma(\vartheta)\simeq\Sigma^{(0)}+\vartheta\Sigma^{(1)}+\vartheta^{2}\Sigma^{(2)}.$
(79)

It is straightforward to note that
$\mathcal{F}(\Sigma(\vartheta),\Sigma(\vartheta))=1$ as a covariance matrix is
always in a full overlap with itself. Another important criteria that is
necessary to impose is,
$\frac{\partial\mathcal{F}(\Sigma(\vartheta),\Sigma(\vartheta+d\vartheta))}{\partial\vartheta}\Bigr{\rvert}_{d\vartheta=0}=0$
[27]. Using the above two conditions, one can expand
$\mathcal{F}(\Sigma(\vartheta),\Sigma(\vartheta+d\vartheta))$ as

$\mathcal{F}(\Sigma(\vartheta),\Sigma(\vartheta+d\vartheta))=1-\frac{\mathcal{F}^{(2)}}{2}d\vartheta^{2}+\mathrm{O}(\vartheta
d\vartheta^{2}+\vartheta^{2}d\vartheta)$ (80)

where $\mathcal{F}^{(2)}=\mathcal{E}^{(2)}+\mathcal{C}^{(2)}$.
$\mathcal{E}^{(2)}$ is proportional to the displacement of the squeezed state
and as a result it is zero. For a single mode scenario, $\mathcal{C}^{(2)}$
has the form

$\begin{split}\mathcal{C}^{(2)}=&\frac{1}{2}\left(\Sigma_{11}^{(0)}\Sigma_{22}^{(2)}+\Sigma_{11}^{(2)}\Sigma_{22}^{(0)}-2\Sigma_{12}^{(0)}\Sigma_{12}^{(2)}\right)\\\
&+\frac{1}{8}\left(\Sigma_{11}^{(1)}\Sigma_{22}^{(1)}-2\Sigma_{12}^{(1)}\Sigma_{12}^{(1)}\right)~{}.\end{split}$
(81)

The quantum Fisher information, in terms of $\mathcal{E}^{(2)}$ and
$\mathcal{C}^{(2)}$, reads [26]

$\mathcal{H}_{\vartheta}=4\mathcal{E}^{(2)}+4\mathcal{C}^{(2)}=4\mathcal{C}^{(2)}$
(82)

as $\mathcal{E}^{(2)}$ is zero for the squeezed bosonic states with no
displacement parameters [14]. In the next part of subsection (IV.2), we shall
obtain the analytical extension of the quantum Fisher information when
quantum-gravity effects are considered in the analysis.

####  IV.2.2 Noise of gravitons and the stochastic average of the quantum
Fisher information

Here, we shall consider the case of the BEC interacting with gravitons. For
the single mode case, the matrix
$\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\tilde{M}}}}_{11}(\varepsilon)$
now has the form

$\begin{split}\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\tilde{M}}}}_{11}(\tilde{\varepsilon})=\begin{bmatrix}1-\frac{\tilde{\varepsilon}}{2\varepsilon}\left[\delta\hat{N}(\mathcal{t})+\frac{\varepsilon\omega_{\beta}\tau\sqrt{\pi}}{4}\left[e^{-\frac{\tau^{2}}{4}(\Omega-2\omega_{\beta})^{2}}-e^{-\frac{\tau^{2}}{4}(\Omega+2\omega_{\beta})^{2}}\right]\right]&&-\frac{\tilde{\varepsilon}}{2\varepsilon}\omega_{\beta}t\delta\hat{N}(\mathcal{t})\\\
\frac{\tilde{\varepsilon}}{2\varepsilon}\omega_{\beta}t\delta\hat{N}(\mathcal{t})&&1+\frac{\tilde{\varepsilon}\omega_{\beta}\tau\sqrt{\pi}}{4}\left[e^{-\frac{\tau^{2}}{4}(\Omega-2\omega_{\beta})^{2}}-e^{-\frac{\tau^{2}}{4}(\Omega+2\omega_{\beta})^{2}}\right]\end{bmatrix}.\end{split}$
(83)

Because of
$\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\mathcal{M}}}}_{11}(\tilde{\varepsilon})$
having elemnts consisting of operators, $\tilde{\Sigma}(\varepsilon)$ from
eq.(68) will also have operator as its elements inspite of
$\Sigma_{\text{sq.}}[0]$ having numbers as its elements. As a result
$4\mathcal{C}^{(2)}$ will now have operatorial contributions in it. As a
result $\mathcal{H}_{\varepsilon}$ will be operator as well. Making use of
eq.(83) in eq.(68), one can obtain the analytical form of
$4\hat{\mathcal{C}}^{2}$ as

$\begin{split}\hat{\mathcal{H}}_{\varepsilon}=&4\hat{\mathcal{C}}^{(2)}=\frac{1}{64}\pi\omega_{\beta}^{2}\tau^{2}\left(e^{2\omega_{\beta}\Omega\tau^{2}}-1\right)^{2}e^{-\frac{\tau^{2}}{2}(\Omega+2\omega_{\beta})^{2}}\left(1+\cosh
4r+(1-3\cos 2\varphi)\sinh^{2}2r\right)\\\
&+\frac{\delta\hat{N}(\mathcal{t})}{32\varepsilon}\left(\sqrt{\pi}\omega_{\beta}\tau\right)\left(e^{2\omega_{\beta}\Omega\tau^{2}}-1\right)e^{-\frac{\tau^{2}}{4}(\Omega+2\omega_{\beta})^{2}}\left(2\cosh^{2}2r+(1-3\cos
2\varphi)\sinh^{2}2r+6\omega_{\beta}\mathcal{t}\sin\varphi\sinh 4r\right)\\\
&+\frac{(\delta\hat{N}(\mathcal{t}))^{2}}{16\varepsilon^{2}}\Bigr{[}1+\cosh
4r-\frac{\sinh^{2}2r}{2}(3+\cos
2\varphi)+\omega_{\beta}\mathcal{t}\bigr{(}\sin\varphi\left(2\sinh^{2}2r(\cos\varphi+\omega_{\beta}\mathcal{t}\sin\varphi)+3\sinh
4r\right)\\\ &+4\omega_{\beta}\mathcal{t}\cosh^{2}2r\bigr{)}\Bigr{]}\\\
=&\mathcal{H}_{\varepsilon}^{(0)}+\frac{\delta\hat{N}(\mathcal{t})}{32\varepsilon}\mathcal{H}_{\varepsilon}^{(1)}+\frac{(\delta\hat{N}(\mathcal{t}))^{2}}{16\varepsilon^{2}}\mathcal{H}_{\varepsilon}^{(2)}\end{split}$
(84)

where

$\displaystyle\mathcal{H}_{\varepsilon}^{(0)}$
$\displaystyle=\frac{1}{64}\pi\omega_{\beta}^{2}\tau^{2}\left(e^{2\omega_{\beta}\Omega\tau^{2}}-1\right)^{2}e^{-\frac{\tau^{2}}{2}(\Omega+2\omega_{\beta})^{2}}\left(1+\cosh
4r+(1-3\cos 2\varphi)\sinh^{2}2r\right)~{},$ (85)
$\displaystyle\mathcal{H}_{\varepsilon}^{(1)}$
$\displaystyle=\left(\sqrt{\pi}\omega_{\beta}\tau\right)\left(e^{2\omega_{\beta}\Omega\tau^{2}}-1\right)e^{-\frac{\tau^{2}}{4}(\Omega+2\omega_{\beta})^{2}}\left(2\cosh^{2}2r+(1-3\cos
2\varphi)\sinh^{2}2r+6\omega_{\beta}\mathcal{t}\sin\varphi\sinh 4r\right)~{},$
(86) $\displaystyle\mathcal{H}_{\varepsilon}^{(2)}$ $\displaystyle=1+\cosh
4r-\frac{\sinh^{2}2r}{2}(3+\cos
2\varphi)+\omega_{\beta}\mathcal{t}\bigr{(}\sin\varphi\left(2\sinh^{2}2r(\cos\varphi+\omega_{\beta}\mathcal{t}\sin\varphi)+3\sinh
4r\right)+4\omega_{\beta}\mathcal{t}\cosh^{2}2r\bigr{)}$ (87)

where $\mathcal{H}_{\varepsilon}^{(0)}$ is the quantum Fisher information for
the classical contribution of the gravitational wave and is exactly similar to
the result obtained in [14]. The other two terms determine the quantum
gravitational contribution to the quantum Fisher information. The quantum
Fisher information operator (or the graviton-noise induced Fisher information)
is not entirely a measurable quantity now. Instead of calling it a Fisher
information operator, it is better to call it a quantum gravitational Fisher
information (QGFI). The straightforward way is to take a stochastic average of
the quantity with respect to the graviton state as

$\begin{split}\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\hat{\mathcal{H}}_{\varepsilon}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}=\mathcal{H}_{\varepsilon}^{(0)}+\frac{\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\\{\delta\hat{N}(\mathcal{t}),\delta\hat{N}(\mathcal{t})\\}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}}{32\varepsilon^{2}}\mathcal{H}_{\varepsilon}^{(2)}\end{split}$
(88)

where the second term from eq.(84) vanishes because the one point correlator
of the noise operator,
$\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\delta\hat{N}(t)\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}$
vanishes. Our next aim is to obtain the analytical form of
$\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\\{\delta\hat{N}(\mathcal{t}),\delta\hat{N}(\mathcal{t})\\}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}$
for the gravitons initially being in a squeezed state. The two point noise
correlator for an arbitrary state of the graviton reads

$\begin{split}&\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\\{\delta\hat{N}_{ij}(t),\delta\hat{N}_{lm}(t^{\prime})\\}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}\\\
&=\frac{4\kappa^{2}}{V}\sum\limits_{\mathbf{k},\mathbf{k}^{\prime}}\sum\limits_{s,s^{\prime}}\epsilon^{s}_{ij}(\mathbf{k})\epsilon^{s^{\prime}}_{lm}(\mathbf{k}^{\prime})\times\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\\{\delta\hat{h}^{s}_{I}(\mathbf{k},t),\delta\hat{h}^{s^{\prime}}_{I}(\mathbf{k}^{\prime},t^{\prime})\\}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}\\\
&=\frac{2\kappa^{2}}{5\pi^{2}}(\delta_{il}\delta_{jm}+\delta_{im}\delta_{jl}-\frac{2}{3}\delta_{ij}\delta_{lm})\smash{\int_{0}^{\Omega_{m}}}dkk^{2}\mathcal{Q}_{\delta
h}(t,t^{\prime},\mathbf{k})\end{split}$ (89)

where in order to obtain the final line of the above equation, we have made
use of the identification between the summation over all possible modes to an
integral over a continuous variable and another definition is used
$\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\\{\delta\hat{h}^{s}_{I}(\mathbf{k},t),\delta\hat{h}^{s^{\prime}}_{I}(\mathbf{k}^{\prime},t^{\prime})\\}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}=\delta_{ss^{\prime}}\delta_{\mathbf{k}+\mathbf{k}^{\prime},0}\mathcal{Q}_{\delta
h}(t,t^{\prime},\mathbf{k})$ [19]. For the graviton initially being in a
squeezed state with squeezing parameter
$r_{k}^{\text{sq.}}=r_{k}e^{i\phi_{k}}$, $\mathcal{Q}_{\delta
h}(t,t^{\prime},\mathbf{k})$ takes the form222For a detailed discusssion on
the graviton state with squeezing, see Appendix (A).

$\begin{split}\mathcal{Q}_{\delta
h}(t,t^{\prime},\mathbf{k})&=\frac{1}{k}\left(\cos\left(k(t-t^{\prime})\right)\cosh
2r_{k}\right.\\\ &\left.-\cos\left(k(t+t^{\prime})-\phi_{k}\right)\sinh
2r_{k}\right)~{}.\end{split}$ (90)

Using eq.(90), one can obtain the two point correlator for $t=t^{\prime}$ as

$\begin{split}&\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\\{\delta\hat{N}_{ij}(t),\delta\hat{N}_{lm}(t)\\}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}=\frac{\kappa^{2}\Omega_{m}^{2}}{5\pi^{2}}\Bigr{(}\delta_{il}\delta_{jm}+\delta_{im}\delta_{jl}\\\
&-\frac{2}{3}\delta_{ij}\delta_{lm}\Bigr{)}\biggr{(}\cosh
2r_{k}+\frac{1}{2\Omega_{m}^{2}t^{2}}\sinh 2r_{k}\bigr{(}\cos\phi_{k}\\\
&-\cos(2\Omega_{m}t-\phi_{k})-2\Omega_{m}t\sin(2\Omega_{m}t-\phi_{k})\bigr{)}\biggr{)}~{}.\end{split}$
(91)

The
$\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\\{\delta\hat{N}(\mathcal{t}),\delta\hat{N}(\mathcal{t})\\}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}$
correlator can now be obtained as

$\begin{split}&\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\\{\delta\hat{N}(\mathcal{t}),\delta\hat{N}(\mathcal{t})\\}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}\equiv\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\\{\delta\hat{N}_{11}(\mathcal{t}),\delta\hat{N}_{11}(\mathcal{t})\\}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}\\\
&=\frac{4\kappa^{2}\Omega_{m}^{2}}{15\pi^{2}}\mathcal{B}(r_{k},\phi_{k},\mathcal{t})\end{split}$
(92)

where the time dependent part of the two-point noise correlator reads

$\begin{split}&\mathcal{B}(r_{k},\phi_{k},\mathcal{t})=\cosh
2r_{k}+\frac{1}{2\Omega_{m}^{2}\mathcal{t}^{2}}\sinh
2r_{k}\bigr{(}\cos\phi_{k}\\\
&-\cos(2\Omega_{m}\mathcal{t}-\phi_{k})-2\Omega_{m}\mathcal{t}\sin(2\Omega_{m}\mathcal{t}-\phi_{k})\bigr{)}~{}.\end{split}$
(93)

The squeezing in the initial graviton states can be very high, even of the
order of $r_{k}\sim 21$ for primordial gravitational wave generated during the
inflationary period [19]. For $\mathcal{t}\rightarrow 0$ limit
$\mathcal{B}(r_{k},\phi_{k},\mathcal{t})$ has the value
$\lim\limits_{\mathcal{t}\rightarrow
0}\mathcal{B}(r_{k},\phi_{k},\mathcal{t})=\cosh 2r_{k}-\cos\phi_{k}\sinh
2r_{k}$ which never vanishes irrespective of any values of $\phi_{k}$. It is
although important to note that for $\phi_{k}=\pi$,
$\mathcal{B}(r_{k},\phi_{k},0)$ becomes maximum and for $\phi_{k}=0$, it
becomes minimum. For a grand unified theory gravitational wave, the cut-off
frequency is at around $\Omega_{m}\sim 10^{8}\text{ Hz}$. In order to observe
the nature of the function $\mathcal{B}$ with respect to time, we use
$\Omega_{m}\sim 10^{8}\text{Hz}$, $r_{k}\sim 10$, and
$\phi_{k}=\\{\frac{\pi}{4},\frac{\pi}{2},\pi\\}$ and plot
$\mathcal{B}(r_{k},\phi_{k},t)$ against $t$ for the above values in Fig.(2).
It is important to note that irrespective of the $\phi_{k}$, in the
$t\rightarrow\infty$ limit it always fluctuates about the same value
$\mathcal{B}(10,\phi_{k},\infty)=\cosh 20\simeq 2.426\times 10^{8}$. It is
easy to observe from Fig.(2) as well as the infinite-time limit that for
$\phi_{k}=\frac{\pi}{2}$,
$\mathcal{B}(r_{k},\frac{\pi}{2},0)=\mathcal{B}(r_{k},\frac{\pi}{2},\infty)$.

Figure 2: We plot for $\phi_{k}=\frac{\pi}{4},\frac{\pi}{2},\pi$ and also plot
the $t\rightarrow\infty$ limit of $\mathcal{B}(r_{k},\phi_{k},t)$ when
$r_{k}=10$ and $\Omega_{m}=10^{8}\text{ Hz}$.

The stochastic average of the quantum Fisher information from eq.(88), after
using eq.(92) takes the form (with proper dimensional reconstruction)

$\begin{split}\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\hat{\mathcal{H}}_{\varepsilon}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}=\mathcal{H}_{\varepsilon}^{(0)}+\frac{\hbar
G\Omega_{m}^{2}}{15\pi\varepsilon^{2}c^{5}}\mathcal{B}(r_{k},\phi_{k},\mathcal{t})\mathcal{H}_{\varepsilon}^{(2)}\end{split}$
(94)

where $\kappa^{2}=\frac{8\pi\hbar G}{c^{3}}$. It is now reasonable to replace
$\mathcal{t}$ by $\tau$ in eq.(94) as the total observation time will be equal
to the single mode measurement time by the BEC. Under this condition, we can
recast eq.(94) as (for $\varphi=\frac{\pi}{2}$)

$\begin{split}\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\hat{\mathcal{H}}_{\varepsilon}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}&=\mathcal{H}_{\varepsilon}^{(0)}+\frac{l_{p}^{2}\Omega_{m}^{2}}{15\pi\varepsilon^{2}c^{2}}\mathcal{B}(r_{k},\phi_{k},\tau)\mathcal{H}_{\varepsilon}^{(2)}\\\
&=\frac{1}{64}\pi\omega_{\beta}^{2}\tau^{2}\left(e^{2\omega_{\beta}\Omega\tau^{2}}-1\right)^{2}e^{-\frac{\tau^{2}}{2}(\Omega+2\omega_{\beta})^{2}}\left(1+\cosh
4r+4\sinh^{2}2r\right)+\frac{l_{p}^{2}\Omega_{m}^{2}}{30\pi\varepsilon^{2}c^{2}}\bigr{(}3+2\omega_{\beta}^{2}\tau^{2}+\cosh
4r+6\omega_{\beta}\tau\\\ &\times\sinh 4r+6\omega_{\beta}^{2}\tau^{2}\cosh
4r\bigr{)}\left(\cosh 2r_{k}+\frac{1}{2\Omega_{m}^{2}\tau^{2}}\sinh
2r_{k}\left(\cos\phi_{k}-\cos(2\Omega_{m}\tau-\phi_{k})-2\Omega_{m}\tau\sin(2\Omega_{m}\tau-\phi_{k})\right)\right)~{}.\end{split}$
(95)

Eq.(95) is one of the main results in our paper. Setting the squeezing angle
to a certain value (here, $\varphi=\frac{\pi}{2}$) is possible and has been
experimentally done [29, 30]. Using eq.(73), we can write down the inequality
in a quantum gravitational setup as

$\begin{split}\langle(\Delta\tilde{\varepsilon})^{2}\rangle\geq\frac{1}{\mathfrak{N}\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\hat{\mathcal{H}}_{\varepsilon}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}}~{}.\end{split}$
(96)

Here, we are considering single mode Bose-Einstein condensate only. As a
result, we can still write down the following relation

$\begin{split}\langle(\Delta\tilde{\varepsilon}_{k_{\beta}})^{2}\rangle=\left(\frac{k_{\beta_{x}}^{2}-k_{\beta_{y}}^{2}}{k_{\beta}^{2}}\right)^{2}\langle(\Delta\varepsilon_{k_{\beta}})^{2}\rangle~{}.\end{split}$
(97)

It is now possible to express $k_{\beta_{x}}$ and $k_{\beta_{y}}$ in the
spherical polar coordinates as
$k_{\beta_{x}}=k_{\beta}\sin\theta_{\beta}\cos\phi_{\beta}$, and
$k_{\beta_{y}}=k_{\beta}\sin\theta_{\beta}\sin\phi_{\beta}$. Using the
spherical polar representation, we obtain
$\langle(\Delta\tilde{\varepsilon}_{k_{\beta}})^{2}\rangle=\sin^{4}\theta_{\beta}\cos^{2}2\phi_{\beta}\langle(\Delta\varepsilon_{k_{\beta}})^{2}\rangle~{}.$
Doing an integral over the first quadrant of the spherical coordinate basis
for all single mode bosonic state of the BEC with momentum $k_{\beta}$, we
obtain $\int
d\Omega_{\beta}=\int_{0}^{\frac{\pi}{2}}d\theta_{\beta}\sin^{5}\theta_{\beta}\int_{0}^{\frac{\pi}{2}}d\phi_{\beta}\cos^{2}2\phi_{\beta}=\frac{2\pi}{15}$.
In particular it is always possible to construct a Bose-Einstein condensate
such that the ground state of the bosonic system consists of only single mode
bosons which is $k_{\beta}$ in the current case. Eq.(96) can then be recast in
the following form

$\langle(\Delta\varepsilon_{k_{\beta}})^{2}\rangle\geq\frac{15}{2\pi\mathfrak{N}\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\hat{\mathcal{H}}_{\varepsilon}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}}.$
(98)

If the time taken for $\mathfrak{N}$ multiple measuerments of the BEC state is
$\mathfrak{t}$ then $\mathfrak{t}\simeq\mathfrak{N}\tau$. At first, we
consider a single measurement of the BEC state which indicates
$\mathfrak{N}=1$. From eq.(98) it is evident that
$\sqrt{\langle(\Delta\varepsilon_{k_{\beta}})^{2}\rangle}$ has the minimum
value at
$\sqrt{\langle(\Delta\varepsilon_{k_{\beta}})^{2}\rangle}_{\text{min.}}=\sqrt{\frac{15}{2\pi\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\hat{\mathcal{H}}_{\varepsilon}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}}}$.
We shall now plot the minimum value of the standard deviation in the amplitude
$\varepsilon$ for a single phonon mode of the BEC against the observation time
$\tau$ in Fig.(3). In order to plot the parameter values used are
$r=0.82,\varphi=\frac{\pi}{2},r_{k}=42,\phi_{k}=\frac{\pi}{2},\Omega=100\text{
Hz}$, and $\omega_{\beta}=50\text{ Hz}$.

Figure 3:
$\sqrt{\langle(\Delta\varepsilon_{k_{\beta}})^{2}\rangle}_{\text{min.}}$ vs
$\tau$ plot for the initial state of the graviton being a highly squeezed
state.

From Fig.(3), it is straightforward to observe that the minimum standard
deviation in the measurement of $\varepsilon_{k_{\beta}}$ corresponding to a
single mode of the BEC is not very high indicating a finite chance of
observation of the graviton. It is although very important to note that with a
decrease in the squeezing of the graviton state, the minimum value in the
measurement of the standard deviation of the gravitational wave amplitude
becomes very high indicating a non detectability of such a scenario. We can
look for the long time behaviour of the
$\sqrt{\langle(\Delta\varepsilon_{k_{\beta}})^{2}\rangle}_{\text{min.}}$ in
Fig.(4) (with same parameters as used to plot Fig.(3)).

Figure 4:
$\sqrt{\langle(\Delta\varepsilon_{k_{\beta}})^{2}\rangle}_{\text{min.}}$ vs
$\tau$ plot to observe the long time behaviour of the minimum value of the
standard deviation of $\varepsilon_{k_{\beta}}$.

We can observe from Fig.(4) that
$\sqrt{\langle(\Delta\varepsilon_{k_{\beta}})^{2}\rangle}_{\text{min.}}$
decreases with increased time of the single measurement of the gravitational
wave indicating a higher chance in detection of the gravity wave. It is
important to note that the graviton signature precisely lies in the detection
of resonance pulses in the single mode BEC initially. We can also check this
analytically. It is important to note that in the $l_{p}\rightarrow 0$ limit,
the contribution from the linearized quantum gravity theory vanishes in
eq.(95) reducing it to the result obtained in [14] which is expected. In this
$l_{p}\rightarrow 0$ limit, if $\tau$ is set to zero, the inequality in
eq.(98) becomes $(\Delta\varepsilon_{k_{\beta}})^{2}\rangle\geq\infty$
implying that no gravitational wave will be detected in such a semiclassical
scenario. The result for the quantum gravitational perspective becomes highly
bizzarre. We observe that

$\lim\limits_{\tau\rightarrow
0}\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\hat{\mathcal{H}}_{\varepsilon}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}=\frac{l_{p}^{2}\Omega_{m}^{2}}{15\pi\varepsilon^{2}c^{2}}\cosh^{2}2r(\cosh
2r_{k}-\cos\phi_{k}\sinh 2r_{k}).$ (99)

For, no squeezing of the initial gravitational wave state ($r_{k}=0$), we
obtain $\lim\limits_{\tau\rightarrow
0}\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\hat{\mathcal{H}}_{\varepsilon}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}=\frac{l_{p}^{2}\Omega_{m}^{2}}{15\pi\varepsilon^{2}c^{2}}\sim
10^{-31}$ ($\Omega_{m}\sim 10^{8}$ and $\varepsilon\sim 10^{-21}$). This
indicates,
$\sqrt{\langle(\Delta\varepsilon)^{2}\rangle}\geq\frac{15\varepsilon
c}{l_{p}\Omega_{m}\sqrt{2}}\sim 10^{16}$. Such a high minimum value of the
$\sqrt{\langle(\Delta\varepsilon)^{2}\rangle}$ parameter indicates a very low
sensitivity of the BEC towards the gravitational wave initially implying an
impossible detection scenario. But things quickly change for a non-vanishing
squeezing of the initial graviton state. Suppose that the initial squeezing
angle is $\phi_{k}=\frac{\pi}{2}$. For
$\langle(\Delta\varepsilon)^{2}\rangle\simeq 1$, the squeezing will be as high
as $r_{k}\simeq 35$. For a grand unified theory inflation $r_{k}\simeq 42$
[19], $\langle(\Delta\varepsilon)^{2}\rangle\simeq 10^{-6}$. This is very
anti-intuitive in a sense that there is a finite possiblity of detecting
primordial gravitational waves from the inflationary time without a proper
time interval of the detector to interact with the gravity wave. In a
linearized quantum gravity model, this is not at all very unphysical as there
is a linearized perturbation field around the BEC even when $\tau=0$. This
will indicate the existence of a gravitons in future generation of BEC based
gravitational wave detection scenario. It is also a possibility that the BEC
itself starts behaving as a gravitating object which may be a very vague
assumption and would not be explored in details in this literature. We leave
this investigation for a future work. Finally, we plot
$\sqrt{\langle(\Delta\varepsilon_{k_{\beta}})^{2}\rangle}_{\text{min.}}$ vs
$\tau$ for various squeezing of the graviton state against the classical case
when no gravitons are present in Fig.(5). For the quantum gravitational case,
we have used the initial squeezing angle of the graviton to be equal to
$\frac{\pi}{2}$. For the BEC state we have used a squeezing of $1.4$ and a
squeezing angle of $\frac{\pi}{2}$ along with the mode frequency is considered
to be at $\omega_{\beta}=50$ Hz. As a result it will be more sensitive for
incoming gravitational wave with frequency $100$ Hz. It is though important to
note that most primordial gravitational waves are supposed to have a frequency
in the $10^{-1}-10$ Hz range.

Figure 5:
$\sqrt{\langle(\Delta\varepsilon_{k_{\beta}})^{2}\rangle}_{\text{min.}}$ vs
$\tau$ plot for squeezed gravitons with squeezing $r_{k}=34$ and $r_{k}=36$
respectively against the case of a classical gravitational wave.

We can observe very important features from Fig.(5). We observe that with the
increase in squeezing
$\sqrt{\langle(\Delta\varepsilon_{k_{\beta}})^{2}\rangle}_{\text{min.}}$
becomes smaller for short measurement periods. It can be seen from the
classical gravitational wave case that
$\sqrt{\langle(\Delta\varepsilon_{k_{\beta}})^{2}\rangle}_{\text{min.}}$
diverges near the $\tau\rightarrow 0$ limit indicating a non-detection of any
gravitational wave. The above result henceforth confirms that for a quantum
gravity scenario, the minimum value of the standard deviation of the gravity
wave amplitude parameter never becomes zero and can be arbitrarily reduced
with squeezed graviton state indicating a higher chance at proving the
existence of gravitons. The next thing that is important to observe is if
there is a standard deviation present in the QGFI. The standard deviation in
the QGFI, reads

$\begin{split}(\Delta\mathcal{H}_{\varepsilon})^{2}&=\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}(\hat{\mathcal{H}}_{\varepsilon}-\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\hat{\mathcal{H}}_{\varepsilon}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}})^{2}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}~{}.\end{split}$
(100)

It is quite straightforward to understand that all odd order correlators will
vanish. Hence, the surving contributions will come from even order
correlators. As the QGFI is calculated upto the second order correlator, we
shall restrict ourselves to second order in the noise correlators only. The
result in eq.(100) reads

$\begin{split}(\Delta\mathcal{H}_{\varepsilon})^{2}&\simeq\frac{\bigr{(}\mathcal{H}_{\varepsilon}^{(1)}\bigr{)}^{2}}{2048\varepsilon^{2}}\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\\{\delta\hat{N}(\tau),\delta\hat{N}(\tau)\\}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}\\\
&\simeq\frac{\hbar
G\Omega_{m}^{2}\bigr{(}\mathcal{H}_{\varepsilon}^{(1)}\bigr{)}^{2}}{960\pi\varepsilon^{2}c^{5}}\mathcal{B}(r_{k},\phi_{k},\tau)~{}.\end{split}$
(101)

The standard deviation in the QGFI can be expressed in an extended form given
as ($\varphi=\frac{\pi}{2}$)

$\begin{split}(\Delta\mathcal{H}_{\varepsilon})^{2}&=\frac{l_{p}^{2}\omega_{\beta}^{2}\Omega_{m}^{2}\tau^{2}}{960\varepsilon^{2}c^{2}}\bigr{(}e^{-\frac{\tau^{2}}{4}(\Omega-2\omega_{\beta})^{2}}-e^{-\frac{\tau^{2}}{4}(\Omega+2\omega_{\beta})^{2}}\bigr{)}^{2}\left(2\cosh^{2}2r+4\sinh^{2}2r+6\omega_{\beta}\tau\sinh
4r\right)^{2}\mathcal{B}(r_{k},\phi_{k},\tau).\end{split}$ (102)

We shall now look at the behaviour of the standard deviation in the QGFI
around resonance. For a finite measurement of $\tau=1$ sec, with
$r_{k}=5,\phi_{k}=\frac{\pi}{2},r=0.83$, and $\varphi=\frac{\pi}{2}$, we plot
$\Delta\mathcal{H}_{\varepsilon}$ against the phonon frequency
$\omega_{\beta}$ for an incoming gravitational wave with frequency 0.1 kHz in
Fig. (6).

Figure 6: $\Delta\mathcal{H}_{\varepsilon}$ vs $\omega_{\beta}$ plot for
$\Omega_{m}=10^{8}$ Hz, $\Omega=100$ Hz, $r_{k}=5$, $\phi_{k}=\frac{\pi}{2}$,
$r=0.83$, and $\varphi=\frac{\pi}{2}$. From the Figure, we can see that the
peak of the standard deviation in the QGFI is observed near the resonance
point $\Omega=2\omega_{\beta}=100$ Hz.

We find out from Fig.(6) that the standard deviation in the stochastic QGFI
will be maximum for the resonance condition which is at
$\Omega=2\omega_{\beta}=100$ Hz. Our next aim is to obtain the QGFI when the
induced noise parameter has a similar decay factor as for the classical
gravitational wave case. What is important to note that the analysis without
any decay factor is much more realistic than the one that we are going to
investigate as any kind of classical decay should not affect the quantum
gravitational influences. At resonance point, with high enough squeezing from
the gravitons, it is possible to enhance the standard deviation of the QGFI to
such an extent that it becomes measurable. Another important aspect can be
obtained by varying the total measurement time-scale $\tau$ for such a
scenario.

Figure 7: $\Delta\mathcal{H}_{\varepsilon}$ vs $\omega_{\beta}$ plot for
$\Omega_{m}=10^{8}$ Hz, $\Omega=100$ Hz, $r_{k}=5$, $\phi_{k}=\frac{\pi}{2}$,
$r=0.83$, and $\varphi=\frac{\pi}{2}$. We have plotted
$\Delta\mathcal{H}_{\varepsilon}$ for different values of $\tau$.

It is very important to observe that the standard deviation in the QGFI
becomes very small with a small measurement time and gets shifted towards the
right for lower values of time $\tau$. This indicates that in order to detect
a standard deviation in the QGFI, one needs to measure it for a longer time
period.

##  V Quantum Fisher information for a decaying noise function

We recall the final solution of the time dependent part of the pseudo-
Goldstone boson in eq.(52) and now use the form of the general noise
fluctuation term at time $t^{\prime}$ as $\delta\hat{N}(t^{\prime})=\cos\Omega
t^{\prime}e^{-\frac{{t^{\prime}}^{2}}{\tau^{2}}}\delta\hat{N}(t).$ The
solution of the time dependent part of the pseudo-Goldstone boson takes the
form for the above noise term as

$\hat{\psi}_{\mathbf{k}_{\beta}}(t)=\hat{\alpha}^{\beta}_{\delta}e^{-i\omega_{\beta}t}+\hat{\beta}^{\beta}_{\delta}e^{i\omega_{\beta}t}$
(103)

where the stochastic Bogoliubov coefficients take the form

$\displaystyle\hat{\alpha}^{\beta}_{\delta}$
$\displaystyle=1-\frac{2\tilde{\varepsilon}\sqrt{\pi}}{\varepsilon}\omega_{\beta}\tau
e^{-\frac{\Omega^{2}\tau^{2}}{4}}\delta\hat{N}(\tau)~{},$ (104)
$\displaystyle\hat{\beta}^{\beta}_{\delta}$
$\displaystyle=\frac{\tilde{\varepsilon}}{4}\sqrt{\pi}\omega_{\beta}\tau\left(e^{-\frac{\tau^{2}}{4}(\Omega-2\omega_{\beta})^{2}}-e^{-\frac{\tau^{2}}{4}(\Omega+2\omega_{\beta})^{2}}\right)$
$\displaystyle+\frac{\tilde{\varepsilon}}{\varepsilon}\sqrt{\pi}\omega_{\beta}\tau\delta\hat{N}(\tau)\left(e^{-\frac{\tau^{2}}{4}(\Omega-2\omega_{\beta})^{2}}-e^{-\frac{\tau^{2}}{4}(\Omega+2\omega_{\beta})^{2}}\right)~{}.$
(105)

One can therefore obtain the stochastic average for the QGFI as follows

$\begin{split}\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\hat{\mathcal{H}}_{\varepsilon}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}=&\mathcal{H}_{\varepsilon}^{(0)}+\frac{\kappa^{2}\omega_{\beta}^{2}\Omega_{m}^{2}\tau^{2}}{15\varepsilon^{2}c^{2}}e^{-\frac{\tau^{2}}{2}(\Omega+2\omega_{\beta})^{2}}\biggr{[}\left[e^{2\omega_{\beta}\Omega\tau^{2}}+1\right]^{2}\\\
\times&(1+\cosh 4r+4\sinh^{2}2r)+4(7-\cosh
4r)e^{2\omega_{\beta}^{2}\tau^{2}}\\\
\times&e^{2\omega_{\beta}\Omega\tau^{2}}\biggr{]}\mathcal{B}(r_{k},\phi_{k},\tau).\end{split}$
(106)

Unlike the previous case, the stochastic average of the QGFI vanishes
initially which is solely due to the decaying nature of the noise parameter.
We shall now sum over all possible modes of the BEC. Before proceeding
further, it is important to note that the BEC was initially quantized within a
box of size $V_{\beta}=L_{\beta}^{3}$ and as a result $k_{\beta}=\frac{\pi
n_{\beta}}{L_{\beta}}$. Eq.(97) can now be recast into the following form

$\begin{split}\langle(\Delta\tilde{\varepsilon})^{2}\rangle&=\sum_{\mathbf{k}_{\beta}}\left(\frac{k_{\beta_{x}}^{2}-k_{\beta_{y}}^{2}}{k_{\beta}^{2}}\right)^{2}\langle(\Delta\varepsilon)^{2}\rangle\\\
\implies\frac{1}{\langle(\Delta\varepsilon)^{2}\rangle}&=\sum_{\mathbf{k}_{\beta}}\left(\frac{k_{\beta_{x}}^{2}-k_{\beta_{y}}^{2}}{k_{\beta}^{2}}\right)^{2}\mathfrak{N}\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\hat{\mathcal{H}}_{\varepsilon}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}\end{split}$
(107)

where in the last line of the above equation, we have made use of the equality
condition from eq.(96) in a quantum gravitational setup. Converting the above
sum over all possible modes to an integration over all possible modes and
defining
$\mathfrak{H}(n_{\beta})\equiv\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\hat{\mathcal{H}}_{\varepsilon}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}$,
we obtain

$\begin{split}\frac{1}{\langle(\Delta\varepsilon)^{2}\rangle}=\frac{2\pi}{15}\int_{0}^{\infty}dn_{\beta}~{}n_{\beta}^{2}\mathfrak{H}(n_{\beta})~{}.\end{split}$
(108)

The inequality in this case can be written as

$\begin{split}\frac{1}{\mathfrak{N}\langle(\Delta\varepsilon)^{2}\rangle}&\leq\frac{c_{s}^{2}\pi^{4}\tau^{2}}{480L_{\beta}^{2}}\mathfrak{r}_{1}\int_{0}^{\infty}dn_{\beta}~{}n_{\beta}^{4}e^{-\frac{(2\pi
c_{s}n_{\beta}+\Omega
L_{\beta})^{2}\tau^{2}}{2L_{\beta}^{2}}}\left(e^{\frac{2\pi
c_{s}n_{\beta}\Omega\tau^{2}}{L_{\beta}}}-1\right)^{2}+\frac{2\kappa^{2}\pi^{2}c_{s}^{2}\Omega_{m}^{2}\tau^{2}}{15\varepsilon^{2}L_{\beta}^{2}}\mathcal{B}(r_{k},\phi_{k},\tau)\\\
&\times\int_{0}^{\infty}dn_{\beta}~{}n_{\beta}^{4}e^{-\frac{(2\pi
c_{s}n_{\beta}+\Omega
L_{\beta})^{2}\tau^{2}}{2L_{\beta}^{2}}}\biggr{[}\mathfrak{r}_{1}\Bigr{(}e^{\frac{2\pi
c_{s}n_{\beta}\Omega\tau^{2}}{L_{\beta}}}+1\Bigr{)}^{2}+\mathfrak{r}_{2}e^{\frac{2\pi
c_{s}n_{\beta}\tau^{2}\left(\Omega+\frac{\pi
c_{s}n_{\beta}}{L_{\beta}}\right)}{L_{\beta}}}\biggr{]}\\\
&=\frac{c_{s}^{2}\pi^{4}\tau^{2}}{480L_{\beta}^{2}}\mathfrak{r}_{1}\mathfrak{I}_{1}+\frac{2\kappa^{2}\pi^{2}c_{s}^{2}\Omega_{m}^{2}\tau^{2}}{15\varepsilon^{2}L_{\beta}^{2}}\mathcal{B}(r_{k},\phi_{k},\tau)\mathfrak{I}_{2}\end{split}$
(109)

where

$\mathfrak{r}_{1}=1+\cosh
4r+4\sinh^{2}2r,~{}\mathfrak{r}_{2}=\frac{1}{2}(7-\cosh 4r)$ (110)

with $\mathfrak{I}_{1}$ and $\mathfrak{I}_{2}$ denoting the first and second
integral in eq.(109). We shall at first explicitly investigate
$\mathfrak{I}_{2}$ as

$\begin{split}\mathfrak{I}_{2}=&\mathfrak{r}_{1}\int_{0}^{\infty}dn_{\beta}~{}n_{\beta}^{4}e^{-\tau^{2}\left[\frac{2\pi^{2}c_{s}^{2}n_{\beta}^{2}}{L_{\beta}^{2}}+\frac{\Omega^{2}}{2}\right]}\cosh^{2}\left[\frac{\pi
c_{s}n_{\beta}\Omega\tau^{2}}{L_{\beta}}\right]\\\
+&\mathfrak{r}_{2}e^{-\frac{\Omega^{2}\tau^{2}}{2}}\int_{0}^{\infty}dn_{\beta}~{}n_{\beta}^{4}~{}.\end{split}$
(111)

The second part of the above integral is divergent. Hence, a way out is to set
the squeezing $r$ of the phonons at such a value that $\mathfrak{r}_{2}$
vanishes effectively. It is straight forward to check that for
$r=\frac{1}{4}\cosh^{-1}(7)\simeq 0.66$, $\mathfrak{r}_{2}$ vanishes. It is
possible to control the squeezing of the phonons. Squeezing of phonons as high
as $r\simeq 0.83$ ($7.2$ dB [31]) has already been achieved and as a result
using a $0.66$ squeezing is of no problem. It is crucial to remember that the
squeezing is generally represented in decibels and the position squeezing $s$
is related to the dimensionless squeezing parameter $r$ via the relation
$s=-10\log_{10}[e^{-2r}]$[32]. In a crystal lattice, using a second order
Raman scattering, phonons have been squeezed [33, 34]. For cold bosonic atoms
in optical lattices [35], such second order Raman scattering [33] or pump-
probe detection scheme [36] can be used to squeeze the phonon modes when an
optical lattice potential is present. One can now obtain obtain the final form
of eq.(109) as

$\begin{split}\frac{1}{\mathfrak{N}\langle(\Delta\varepsilon)^{2}\rangle}\leq\frac{V_{\beta}\mathfrak{r}_{1}}{7680\sqrt{2\pi}c_{s}^{3}\tau^{3}}(\Omega^{4}\tau^{4}+6\Omega^{2}\tau^{2}+3-3e^{-\frac{\Omega^{2}\tau^{2}}{2}})+\frac{\hbar
GV_{\beta}\Omega_{m}^{2}\mathfrak{r}_{1}\mathcal{B}(\tau)}{225\pi\sqrt{2\pi}\varepsilon^{2}c_{s}^{3}\tau^{3}c^{5}}(\Omega^{4}\tau^{4}+6\Omega^{2}\tau^{2}+3+3e^{-\frac{\Omega^{2}\tau^{2}}{2}})\end{split}$
(112)

where we have used $\mathcal{B}(\tau)\equiv\mathcal{B}(r_{k},\phi_{k},\tau)$.
For the right hand side of the above equation no apprximation for the
$\Omega\tau$ factor has been taken. Now, we shall investigate into the case
when $\Omega\sim\Omega_{m}$ ($10^{8}$ Hz). It is important to note that the
BEC in general is prepared in a single length direction and the perpendicular
directions are quite smaller. Our model on the other hand carries the cubic
BEC approximation. It has been possible to create a BEC with length
$L_{\beta}\sim 10^{-3}$ m [38, 39, 37]. As $\tau$ is the duration of the
single measurement of the gravitational wave
$\tau=\frac{v_{\text{max}}}{L_{\beta}}=\frac{c}{L_{\beta}}\sim 10^{-11}$ sec.
Hence, for $\Omega\sim\Omega_{m}$, $\Omega\tau\sim 10^{-3}$. As a result,
$\Omega\tau\ll 1$. For a total observation time of $\tau_{\text{obs.}}$, one
can run approximately $\mathfrak{N}\sim\frac{\tau_{\text{obs.}}}{\tau}$ number
of observations. Under the $\Omega\tau\ll 1$ condition, eq.(112) can be recast
as

$\begin{split}\langle(\Delta\varepsilon)^{2}\rangle\gtrsim&\frac{1024\sqrt{2\pi}c_{s}^{3}\tau^{2}}{\Omega^{2}V_{\beta}\tau_{\text{obs.}}\mathfrak{r}_{1}}\biggr{(}1-\frac{1024l_{p}^{2}\Omega_{m}^{2}}{50\pi\varepsilon^{2}c^{2}}\mathcal{B}(\tau)\\\
&-\frac{2048l_{p}^{2}\Omega_{m}^{2}}{75\pi\varepsilon^{2}c^{2}\Omega^{2}\tau^{2}}\mathcal{B}(\tau)\biggr{)}~{}.\end{split}$
(113)

It is important to note that $(\Delta\varepsilon)^{2}$ can never be negative,
as a result from the equality condition we can write down the minimum value
for the observation time of the single measurement of the gravitational wave
$\tau$ to be

$\begin{split}\tau_{\text{min.}}\simeq\sqrt{\frac{2}{3\pi}}\frac{32l_{p}\Omega_{m}}{c\varepsilon\Omega}\mathcal{B}(\tau)~{}.\end{split}$
(114)

For a vacuum state without any squeezing and $\Omega\sim\Omega_{m}$, $\tau$
attains its absolute minimum value, which is given by

$\begin{split}\tau_{\text{min.}}^{0}\simeq&\sqrt{\frac{2}{3\pi}}\frac{32l_{p}\Omega_{m}}{c\varepsilon\Omega}\Bigr{\rvert}_{\Omega\rightarrow\Omega_{m}}\sim\frac{1.59\times
10^{-14}}{\Omega_{m}}\text{ sec}\\\
\implies\tau_{\text{min.}}^{0}\simeq&1.59\times 10^{-22}\text{
sec}.\end{split}$ (115)

This is a very important result in our paper. In [14], it was argued that the
measurement cannot be arbitrarily smaller by comparing numerical data. In our
case, a complete quantum gravity calculation puts up a theoretical lower bound
for the measurement time when the noise fluctuation is weighted by a Gaussian
decay factor. Eq.(115) reveals that the single measurement time $\tau$ must be
greater than or equal to $\tau^{0}_{\text{min.}}$. It is although very
important to note that the Gaussian decay term in the classical part of the
gravitational waves comes entirely from the template of the gravitational wave
whereas in this section it is imposed by hand. Therefore, the results obtained
in section (IV) are much more plausible than this one. Although one indeed can
induce such Gaussian decay mechanically into the system which shall lead to a
much more complicated result than the simpler model presented here.

###  V.1 BEC as a graviton detector

In this subsection, we shall argue that the BEC will suffice as a graviton
detector when future generation of gravitational wave detector will come up.
We shall here use the projected sensitivity of the upcoming LISA333The full
form of LISA is Laser Interferometer Space Antenna. observatory as a baseline
for the comparison. In the next section, we shall consider a more realistic
case when there are interaction between the phonon modes which will result in
a decoherence effect. We start with the sensitivity formula presented in
Science Requirement Document (SciRD) [40] projected for the LISA observatory.
A detailed discussion can be obtained in [41]. The SciRD sensitivity formula
reads [40, 41]

$\displaystyle S_{\text{h,SciRD}}(f)$
$\displaystyle=\frac{10}{3}\left(\frac{S_{I}(f)}{(2\pi
f)^{4}}+S_{II}(f)\right)R(f)\text{Hz}^{-1}$ (116) $\displaystyle S_{I}(f)$
$\displaystyle=5.76\times
10^{-48}\left(1+\frac{f_{1}^{2}}{f^{2}}\right)\text{sec}^{-4}\cdot\text{Hz}^{-1}$
(117) $\displaystyle S_{II}(f)$ $\displaystyle=3.6\times
10^{-41}\text{Hz}^{-1}$ (118) $\displaystyle R(f)$
$\displaystyle=1+\frac{f^{2}}{f_{2}^{2}}$ (119)

where $f_{1}=0.4\text{mHz}$ and $f_{2}=25\text{mHz}$. In order to compare the
above result we consider the equality from eq.(112) and use the minimum value
of the standard deviation in the amplitude parameter
$\sqrt{\langle(\Delta\varepsilon)^{2}\rangle}_{\text{min}}\bigr{\rvert}_{\Omega=f}$.
The sensitivity of the BEC is given by
$\frac{\sqrt{\langle(\Delta\varepsilon)^{2}\rangle}_{\text{min}}}{\sqrt{f}}$
$\text{Hz}^{-\frac{1}{2}}$. We use the following parameter values
$\tau=10^{-6}\text{ sec}$, $\tau_{\text{obs}}=10^{2}\text{ sec}$,
$L_{\beta}=10^{-3}\text{ m}$, and $c_{s}=0.012\text{ m}\cdot\text{sec}^{-1}$.
The plot of
$\frac{\sqrt{\langle(\Delta\varepsilon)^{2}\rangle}_{\text{min}}}{\sqrt{f}}$
for the BEC vs the SciRD sensitivity formula
($\sqrt{S_{\text{h,SciRD}}(f)}\text{ Hz}^{-\frac{1}{2}}$) is plotted in Fig.
(8). It is important to note that LISA is mainly going to work for the
detection of very low frequency gravitational waves (especially primordial
gravitational waves).

Figure 8: The SciRD sensitivity formula is plotted along with the BEC-graviton
model sensitivity formula against the wave frequency $f$. The BEC sensitivity
plot is viable when the resonance condition is satisfied which is
$2\omega_{\beta}=\Omega=f$.

From Fig.(8), it is evident that with higher squeezing from the gravitons
lower frequencies can be probed by the BEC. It is important to note that SciRD
plot for the LISA targets classical gravitational waves. Hence, a simultaneous
detection by LISA and a BEC will prove the existence of gravitons. We now
compare this SciRD sensitivity formula with the case when the BEC is
interacting as a classical gravitational wave in Fig.(9).

Figure 9: The SciRD sensitivity formula is plotted along with the BEC model
sensitivity formula against the wave frequency $f$. We plot the case of BEC-
Classical Gravity wave model alongside the BEC-Graviton model.

In Fig.(9), we observe that the semi-classical BEC model with classical
gravity wave interaction is not a good candidate for detecting low-frequency
gravitational waves. On the contrary with same phonon squeezing the BEC can
detect graviton signatures when the gravitons are arriving with high enough
squeezing. This reinforces our result, that a BEC is one of the best
candidates for capturing signature of gravitons. When a gravity wave in the
common frequency range will be detected by LISA, a detection by a BEC confirms
the fact that gravitons exist as a classical gravity wave will never be
detected by a BEC in such low-frequency ranges with fixed squeezing as low as
$r=2.3$. It also confirms that a BEC will better serve as a graviton detector
than a classical gravity wave detector.

##  VI Effects of Decoherence from interacting modes of the phonons on the
“QGFI”

Upto the previous section, we have considered dissipative system only. In this
section, interaction among the phonon modes will be considered as a result
there will be dissipation inside of the system. The simple idea is to connect
the single mode BEC system with a thermal bath. For single mode Gaussian
state, the time evolution of the covariance matrix reads

$\begin{split}\Sigma(t)=\Gamma(t)\Sigma_{0}\Gamma^{T}(t)+\Sigma_{\infty}(t)\end{split}$
(120)

where $\Sigma_{0}$ is the covariance matrix of the single mode Gaussian state
initially and $\Gamma(t)=e^{-\frac{\gamma t}{2}}\mathbb{1}_{2}$ with $\gamma$
being the dissipation constant. Here in eq.(120), $\Sigma_{\infty}(t)$ denotes
the time dependent covariance matrix of the Gaussian reservoir and is given by

$\Sigma_{\infty}(t)=(1-e^{-\gamma t})\Sigma_{\infty}.$ (121)

Using eq.(s)(120,121), one can write down the elements of the covariance
matrix as [42, 43]

$\Sigma_{ij}(t)=e^{-\gamma t}{\Sigma_{0}}_{ij}+(1-e^{-\gamma
t}){\Sigma_{\infty}}_{ij}$ (122)

where $i,j=1,2$. In this section, we have followed the analysis presented in
[12, 43]. The purity of the quantum state is given by
$\mu(t)=\frac{1}{2\sqrt{\det[\Sigma(t)]}}$. In such a scenario, the elements
of the covariance matrix takes the form

$\begin{split}\Sigma_{11}(t)&=\frac{1}{2\mu(t)}\left(\cosh
2r(t)+\cos\varphi(t)\sinh 2r(t)\right)~{},\\\
\Sigma_{12}(t)&=\Sigma_{21}(t)=\frac{1}{2\mu(t)}\sin\varphi(t)\sinh
2r(t)~{},\\\ \Sigma_{22}(t)&=\frac{1}{2\mu(t)}\left(\cosh
2r(t)-\cos\varphi(t)\sinh 2r(t)\right)~{}.\end{split}$ (123)

In the above equation the squeezing parameter and squeezing phase both becomes
time dependent due to dissipation in the system. One can use
$\mu(0)\equiv\mu_{0}$ and $r(t)\equiv r_{0}$ to define the purity and
squeezing initially of the single mode bosonic system and the elements of the
initial covariance matrix read

$\begin{split}{\Sigma_{0}}_{11}&=\frac{1}{2\mu_{0}}\left(\cosh
2r_{0}+\cos\varphi_{0}\sinh 2r_{0}\right)~{},\\\
{\Sigma_{0}}_{12}&={\Sigma_{0}}_{21}=\frac{1}{2\mu_{0}}\sin\varphi_{0}\sinh
2r_{0}~{},\\\ {\Sigma_{0}}_{22}&=\frac{1}{2\mu_{0}}\left(\cosh
2r_{0}-\cos\varphi_{0}\sinh 2r_{0}\right)\end{split}$ (124)

The covariance matric elements of the Gaussian reservoir initially reads

$\begin{split}{\Sigma_{\infty}}_{11}&=\frac{1}{2\mu_{\infty}}\left(\cosh
2r_{\infty}+\cos\varphi_{\infty}\sinh 2r_{\infty}\right)~{},\\\
{\Sigma_{\infty}}_{12}&={\Sigma_{\infty}}_{21}=\frac{1}{2\mu_{\infty}}\sin\varphi_{\infty}\sinh
2r_{\infty}~{},\\\ {\Sigma_{\infty}}_{22}&=\frac{1}{2\mu_{\infty}}\left(\cosh
2r_{\infty}-\cos\varphi_{\infty}\sinh 2r_{\infty}\right)\end{split}$ (125)

where $\mu_{\infty}$, $r_{\infty}$ and $\varphi_{\infty}$ denote respectively
the purity, squeezing parameter, and squeezing angle of the reservoir. One can
easily consider a thermal bath with no squeezing which is given by the
condition $r_{\infty}=0$ [12, 42] and reduces the covariance matrix of the
thermal bath to $\Sigma_{\infty}=\frac{1}{2\mu_{\infty}}\mathbb{1}_{2}$.
Initially, we shall start with non-zero squeezing for the thermal bath and
later will reduce down to the no squeezing case. Using eq.(123,124,125) in
eq.(122), one obtains three equations which are given by

$\displaystyle\frac{\cosh 2r(t)+\cos\varphi(t)\sinh 2r(t)}{2\mu(t)}$
$\displaystyle=\frac{e^{-\gamma t}}{2\mu_{0}}\left(\cosh
2r_{0}+\cos\varphi_{0}\sinh 2r_{0}\right)+\frac{1-e^{-\gamma
t}}{2\mu_{\infty}}\left(\cosh 2r_{\infty}+\cos\varphi_{\infty}\sinh
2r_{\infty}\right)$ (126) $\displaystyle\frac{1}{2\mu(t)}\sin\varphi(t)\sinh
2r(t)$ $\displaystyle=\frac{e^{-\gamma t}}{2\mu_{0}}\sin\varphi_{0}\sinh
2r_{0}+\frac{1-e^{-\gamma t}}{2\mu_{\infty}}\sin\varphi_{\infty}\sinh
2r_{\infty}$ (127) $\displaystyle\frac{\cosh 2r(t)-\cos\varphi(t)\sinh
2r(t)}{2\mu(t)}$ $\displaystyle=\frac{e^{-\gamma t}}{2\mu_{0}}\left(\cosh
2r_{0}-\cos\varphi_{0}\sinh 2r_{0}\right)+\frac{1-e^{-\gamma
t}}{2\mu_{\infty}}\left(\cosh 2r_{\infty}-\cos\varphi_{\infty}\sinh
2r_{\infty}\right).$ (128)

Using the above three equations, one obtains three equations describing the
dissipiation relations of the three independent parameters as [43]

$\displaystyle\mu(t)=\mu_{0}\left[e^{-2\gamma
t}+\frac{\mu_{0}^{2}}{\mu_{\infty}^{2}}(1-e^{-\gamma
t})^{2}+\frac{2\mu_{0}}{\mu_{\infty}}e^{-\gamma t}(1-e^{-\gamma t})\left(\cosh
2r_{0}\cosh 2r_{\infty}-\cos(\varphi_{0}-\varphi_{\infty})\sinh 2r_{0}\sinh
2r_{\infty}\right)\right]^{-\frac{1}{2}},$ (129) $\displaystyle\cosh
2r(t)=\mu(t)\left[\frac{e^{-\gamma t}}{\mu_{0}}\cosh 2r_{0}+\frac{1-e^{-\gamma
t}}{\mu_{\infty}}\cosh
2r_{\infty}\right],~{}\tan\varphi(t)=\frac{\sin\varphi_{0}\sin
2r_{0}+\frac{\mu_{0}}{\mu_{\infty}}(e^{\gamma t}-1)\sin\varphi_{\infty}\sinh
2r_{\infty}}{\cos\varphi_{0}\sin 2r_{0}+\frac{\mu_{0}}{\mu_{\infty}}(e^{\gamma
t}-1)\cos\varphi_{\infty}\sinh 2r_{\infty}}~{}.$ (130)

Our calculation produces results slightly different than the one presented in
[43]. One of the primary reasons is that the different signature of the off-
diagonal elements of the covariance matrix corresponding to the single mode
Gaussian state of the Bose-Einstein condensate. We shall now set
$r_{\infty}=0$ which recasts the ‘$\varphi(t)$’ equation in eq.(130) to
$\tan\varphi(t)=\tan\varphi_{0}$. This implies that the squeezing angle does
not change overtime if the reservoir attached has no squeezing. Hence, we can
replace $\phi(t)$ by $\phi_{0}$ in our analysis. Eq.(s)(129,130), in this non-
squeezed thermal bath consideration, then reduces to [14, 42]

$\begin{split}\mu(t)=&\mu_{0}\biggr{(}e^{-2\gamma
t}+\frac{\mu_{0}^{2}}{\mu_{\infty}^{2}}(1-e^{-\gamma
t})^{2}+\frac{2\mu_{0}}{\mu_{\infty}}e^{-\gamma t}(1-e^{-\gamma t})\cosh
2r_{0}\biggr{)}^{-\frac{1}{2}},~{}\cosh 2r(t)=\mu(t)\left[\frac{e^{-\gamma
t}}{\mu_{0}}\cosh 2r_{0}+\frac{1-e^{-\gamma
t}}{\mu_{\infty}}\right]~{}.\end{split}$ (131)

Assuming that
$r_{0}>\max\left[\frac{\mu_{0}}{\mu_{\infty}},\frac{\mu_{\infty}}{\mu_{0}}\right]$,
one can get the value of $t$ for which the purity becomes minimum as [42, 43]

$t_{\text{min.}}=\frac{1}{\gamma}\ln\left[\frac{\mu_{0}^{2}+\mu_{\infty}^{2}-2\mu_{0}\mu_{\infty}\cosh
2r_{0}}{\mu_{0}^{2}-\mu_{0}\mu_{\infty}\cosh 2r_{0}}\right]~{}.$ (132)

Here, $t_{\text{min.}}$ serves as the characteristic decoherence time of the
squeezed single-mode bosonic states. It is straightforward to understand that
the way to incorporate the dissipation into the theory is to replace $r$ by
$r(t)$ into the stochastic average of the QGFI from eq.(95). We can rewrite
the stochastic average of the QGFI from eq.(95) as

$\begin{split}\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\hat{\mathcal{H}}_{\varepsilon}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}&=\frac{1}{32}\pi\omega_{\beta}^{2}\tau^{2}\left(e^{2\omega_{\beta}\Omega\tau^{2}}-1\right)^{2}e^{-\frac{\tau^{2}}{2}(\Omega+2\omega_{\beta})^{2}}\left(3\cosh^{2}2r(\tau)-2\right)+\frac{l_{p}^{2}\Omega_{m}^{2}}{15\pi\varepsilon^{2}c^{2}}\Bigr{(}1-2\omega_{\beta}^{2}\tau^{2}+\cosh^{2}2r(\tau)(1+6\omega_{\beta}^{2}\tau^{2})\\\
&+6\omega_{\beta}\tau\cosh
2r(\tau)\sqrt{\cosh^{2}2r(\tau)-1}\Bigr{)}\mathcal{B}(r_{k},\phi_{k},\tau)\end{split}$
(133)

where we have replaced $t$ by $\tau$ in $r(t)$. The simplest way to
incorporate dissipation into the theory is by replacing the $\cosh 2r(\tau)$
terms using eq.(131). Instead of doing an analytical calculation, we need to
compare the result using plots. It is important to note that Beliaev damping
will be dominant at low temperature which in the zero temperature limit takes
the form [46, 44, 45]

$\gamma\simeq\frac{3}{640\pi}\frac{\hbar\omega_{\beta}^{5}}{m_{\beta}n_{\beta}c_{s}^{5}}$
(134)

where $n_{\beta}$ denotes the number density of the atoms in the BEC and
$m_{\beta}$ denotes the mass of each individual atoms. For a Bose-Einstein
condensate with a number density of $7\times 10^{20}\text{ m}^{-3}$,
$c_{s}\simeq 1.2\times 10^{-2}\text{ m sec}^{-1}$[47].

Figure 10:
$\sqrt{\langle(\Delta\varepsilon_{k_{\beta}})^{2}\rangle}_{\text{min.}}$ vs
$\tau$ has been plotted for the case with and without damping with $r_{k}=33$.
Figure 11: The minimum value for the standard deviation in the amplitude
parameter $\varepsilon_{k_{\beta}}$ has been plotted against the observation
time $\tau$ for the case with and without damping with $r_{k}=2$.

For $\omega_{\beta}=10$ Hz, one obtains $\gamma\simeq 9.034\times
10^{-19}\text{ sec}^{-1}$. We consider the initial state of the BEC as well as
the thermal bath to be pure ($\mu_{0}=\mu_{\infty}=1$). We consider mainly two
cases. The first case when the squeezing of the graviton state is $r_{k}=33$
and the second case $r_{k}=2$. Both of the squeezing angles here, are set to
$\frac{\pi}{2}$, the squeezing if the phonon is at $r=0.83$ and the incoming
gravitational wave has a frequency $\omega=20$ Hz. For the first case we do
not observe any difference due to damping untill a very later time where the
decoherence results in a faster decay of the
$\sqrt{\langle(\Delta\varepsilon_{k_{\beta}})^{2}\rangle}_{\text{min.}}$ with
the observation time as can be seen from Fig.(10). For Fig.(11), the squeezing
is reduced to $r_{k}=2$ and as a result
$\sqrt{\langle(\Delta\varepsilon_{k_{\beta}})^{2}\rangle}_{\text{min.}}$
separates out at a very early observation time $\tau$ ($\sim 7\times 10^{-7}$
sec). This indicates that decoherence effect becomes way less significant for
a lower squeezing of the initial graviton state. One can also investigate the
decoherence effect for the case of the decaying noise function presented in
section (V). For a correct incorporation of decoherence due to interacting
phonon modes in the theory, one should follow the prescription in [48].
Another important point to note is that in order to conduct such metrological
measurements over a time period $\tau_{\text{obs.}}$, one needs to
continuously generate the Bose-Einstein condensate using magneto-optical
traps, the experimental setup of which has been proposed in [49] and later
observed in [50]. It is important to note that we are mainly looking for
signatures of quantum gravity using a Bose-Einstein condensate. In the next
part of this paper we shall delve into the fundamental effects of a linearized
quantum gravity theory on a BEC.

##  VII Conclusion

In this paper, we have considered the simplest model of a Bose-Einstein
condensate interacting with an incoming gravitational wave. The background is
considered to be a flat Minkowski spacetime with fluctuations over it. In
order to incorporate quantum gravitational effects into the theory, we have
quantized the gravitational perturbation, over the flat background by doing a
discrete Fourier mode decomposition, raising the phase space variables to
operator status and applying suitable commutation relation between the
conjugate variables. Using the principle of least action, an equation of
motion both corresponding to the graviton as well as the time dependent part
of the pseudo-Goldstone boson is obtained. Because of the involvement of the
graviton part, the equation of motion corresponding to the boson becomes
stochastic or Langevin-like in nature. As a result the solution obtained for
the bosonic modes become stochastic as well. It is important to note that the
Bogoliubov coefficients obtained in [14] now have contributions from the noise
fluctuation which raises the coefficients to an operator status. We have then
used quantum metrological techniques to obtain the quantum Fisher information.
Because of the quantum gravitational analysis, the quantum Fisher information
picks up effects from the noise fluctuation making it stochastic in nature.
This quantum gravity modified Fisher information is completely a new quantity
and is termed as the “quantum gravitational Fisher information” (QGFI) in this
paper. It is evident that the observable will be the stochastic average of the
QGFI. It is important to note that we have used squeezed graviton states. The
QGFI gives us very fundamental insigts into the detection scenario. From the
Fig.(s)(3,4,5) and the analytical calculation, we observe that with high
enough squeezing there is a finite probabillity of the detection of a graviton
background even for a very small observation time. This is not a very bizarre
scenario as for a quantum gravity consideration, there is a background field
always present. Hence, if a very small measurement using a single mode BEC can
be done just initially, it will definitely be a graviton signature especially
due to squeezed gravitons. This is the main result of our paper and it will
lead completely towards a new era of graviton detection models using BEC. We
have then calculated some other important aspects of such a model. We have
calculated the standard deviation in the QGFI and observe that it maximizes
for a higher observation time at the resonance condition. It is also possible
to measure the the standard deviation of the QGFI. Although it will be very
complicated, we hope for such observations in advanced experimental scenarios
using continuously generated Bose-Einstein condensates. Next, we have
considered a different scenario where the external noise fluctuation gets
attenuated overtime by the use of a Gaussian decay factor. This analysis helps
us to decay out the unusual noise fluctuations created due to the noise of
gravitons. This analysis is very important in a sense that it helps us to
obtain an absolute lower bound to the time of single measurement of the
gravitational wave and is of the order of $10^{-22}$ sec. Next, we have used
the required LISA sensitivity fromula [40, 41] and comparing with our results,
we find that a BEC will be one of the best candidates for a graviton
detection. In order to detect a graviton, the graviton must come with high
enough squeezing which can only exist in primordial gravitational waves coming
from the inflationary time period. This is another very important observation
in our paper which shows that even without high phonon squeezing [14, 15], the
BEC will act as a graviton detector. Finally, we have considered a more
realistic scenario when the phonon modes of the Bose-Einstein condensates are
interacting. We have reproduced the results of the time dependence of the
purity of the Bose-Einstein condensate as well as the phonon-squeezing
parameter for the covariance matrix obtained in our case. Finally, we have
obtained the form of stochastic average of the QGFI when decoherence is
present in the theory. In order to truly observe the behaviour of the minimum
value of the standard deviation in the amplitude $\varepsilon_{k_{\beta}}$, we
have plotted it against the single observation time for the case when
decoherence is present and when decoherence is not present in the system. We
have plotted for the cases of a high graviton squeezing and low graviton
squeezing. It is important to note that the change in the minimum standard
deviation $\varepsilon_{k_{\beta}}$ becomes way less significant at initial
times for high enough squeezing. For almost very small (even for no squeezing
case), the difference becomes significant even at initial times but the
standard deviation value suggests (Fig.(11)) that such effects will not be
observable at such initial times. This is a very complicated experimental
scenario and will be very difficult to perform as the detection of the
graviton signatures realizes highly on the accuracy of instantaneous
measurement. Hence, the way out is to make multiple measurements and if a
resonance spike is observed in the pico-nano second time regime (even micro
second) from the starting of a single measurement, it shall be a conclusive
evidence of the existence of a graviton. In the next part $``\textit{Zweite
Abhandlung}"$, we shall explore a much more fundamental scenario where quantum
gravity will play a leading role and upon experimental verification will be a
conclusive evidence of the quantum nature of gravity (specifically the
evidence of linearized quantum gravity)444This analysis is an extended version
of the letter [51]..

##  Appendix A Squeezed graviton state and the two point correlator

In this appendix, we shall calculate the two point correlator for the initial
graviton state to be in a squeezed state and try to obtain eq.(90) in this
process. The initial graviton is considered to be in a squeezed state. If the
squeezing and the displacement operators are given by

$\displaystyle\hat{S}(r^{\text{sq.}})$
$\displaystyle=e^{\frac{1}{V}\sum\limits_{\mathbf{k},s}\left(r^{\text{sq.}*}_{k}\hat{a}_{s}(\mathbf{k})\hat{a}_{s}(-\mathbf{k})+r^{\text{sq.}}_{\mathbf{k}}\hat{a}_{s}^{\dagger}(\mathbf{k})\hat{a}_{s}^{\dagger}(-\mathbf{k})\right)},$
(135) $\displaystyle\hat{D}(\mathfrak{B})$
$\displaystyle=e^{\frac{1}{V}\sum\limits_{\mathbf{k},s}\left(\mathfrak{B}_{k}\hat{a}^{\dagger}_{s}(\mathbf{k})-\mathfrak{B}_{k}\hat{a}_{s}(\mathbf{k})\right)}$
(136)

where $r^{\text{sq.}}_{k}=r_{k}e^{i\phi_{k}}$, then the displaced squeezed
state reads

$|r^{\text{sq.}},\mathfrak{B}\rangle=\hat{S}(r^{\text{sq.}})\hat{D}(\mathfrak{B})|0\rangle.$
(137)

If the Minkowski mode solution is given by
$u_{k}(t)=\frac{1}{\sqrt{2k}}e^{-ikt}$, then the squeezed mode function has
the form

$u_{k}^{\text{sq.}}(t)=u_{k}(t)\cosh r_{k}-e^{-i\phi_{k}}u_{k}^{*}(t)\sinh
r_{k}.$ (138)

Our primary aim is to calculate
$\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\\{\delta\hat{h}_{I}^{s}(\mathbf{k},t),\delta\hat{h}_{I}^{s^{\prime}}(\mathbf{k}^{\prime},t)\\}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}$
where the expectation is taken with respect to the state in eq.(137). We
already know that
$\hat{h}^{s}_{I}(\mathbf{k},t)=\hat{a}_{s}(\mathbf{k})u_{k}(t)+\hat{a}^{\dagger}_{s}(-\mathbf{k})u_{k}^{*}(t)$
and $\hat{\delta
h}_{I}^{s}(\mathbf{k},t)=\hat{h}_{I}^{s}(\mathbf{k},t)-\langle\hat{h}_{I}^{s}(\mathbf{k},t)\rangle.$
Before proceeding further, we want to write down the following two relations
(note that both the $\hat{S}$ and $\hat{D}$ operators are unitary)

$\displaystyle\hat{\mathcal{A}}_{\mathbf{k},s}(r^{\text{sq.}},\mathfrak{B})=\hat{D}^{\dagger}(\mathfrak{B})\hat{S}^{\dagger}(r^{\text{sq.}})\hat{a}_{s}(\mathbf{k})\hat{S}(r^{\text{sq.}})\hat{D}(\mathfrak{B}).$
(139)

It is then straightforward to obtain the following two relations by making use
of eq.(s)(135,136) as

$\displaystyle\hat{\mathcal{A}}_{\mathbf{k},s}(r^{\text{sq.}},\mathfrak{B})$
$\displaystyle=$ $\displaystyle(\hat{a}_{s}(\mathbf{k})+\mathfrak{B}_{k})\cosh
r_{k}-(\hat{a}^{\dagger}_{s}(-\mathbf{k})+\mathfrak{B}^{*}_{k})e^{i\phi_{k}}\sinh
r_{k}~{},$ (140)
$\displaystyle\hat{\mathcal{A}}^{\dagger}_{-\mathbf{k},s}(r^{\text{sq.}},\mathfrak{B})$
$\displaystyle=$
$\displaystyle(\hat{a}_{s}^{\dagger}(-\mathbf{k})+\mathfrak{B}^{*}_{k})\cosh
r_{k}-(\hat{a}_{s}(\mathbf{k})+\mathfrak{B}_{k})e^{-i\phi_{k}}\sinh r_{k}$
(141)

where we have made use of the fact that the sign of $k$ remains invariant for
both $\mathbf{k}$ and $\mathbf{-k}$. Using the above two relations, we obtain

$\begin{split}&\hat{D}^{\dagger}(\mathfrak{B})\hat{S}^{\dagger}(r^{\text{sq.}})\hat{\delta
h}_{I}^{s}(\mathbf{k},t)\hat{S}(r^{\text{sq.}})\hat{D}(\mathfrak{B})\\\
&=u_{k}(t)\hat{\mathcal{A}}_{\mathbf{k},s}(r^{\text{sq.}},\mathfrak{B})+u_{k}^{*}(t)\hat{\mathcal{A}}^{\dagger}_{-\mathbf{k},s}(r^{\text{sq.}},\mathfrak{B})-u_{k}(t)\langle\hat{\mathcal{A}}_{\mathbf{k},s}(r^{\text{sq.}},\mathfrak{B})\rangle-
u_{k}^{*}(t)\langle\hat{\mathcal{A}}^{\dagger}_{-\mathbf{k},s}(r^{\text{sq.}},\mathfrak{B})\rangle\\\
&=u_{k}^{\text{sq.}}(t)\hat{a}_{s}(\mathbf{k})+u_{k}^{\text{sq.}*}(t)\hat{a}_{s}^{\dagger}(-\mathbf{k})\end{split}$
(142)

where we have made use of eq.(138) to arrive at the final line of the above
equation. We already know the commutation relation among the ladder operators
as
$[\hat{a}_{s}(\mathbf{k}),\hat{a}^{\dagger}_{s}(-\mathbf{k})]=\delta_{s,s^{\prime}}\delta_{\mathbf{k},-\mathbf{k}^{\prime}}$.
Using the above results one obtains the following relation for the two-point
correlator as

$\begin{split}\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\\{\delta\hat{h}_{I}^{s}(\mathbf{k},t),\delta\hat{h}_{I}^{s^{\prime}}(\mathbf{k}^{\prime},t)\\}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}&=\langle
r^{\text{sq.}},\mathfrak{B}|\\{\delta\hat{h}_{I}^{s}(\mathbf{k},t),\delta\hat{h}_{I}^{s^{\prime}}(\mathbf{k}^{\prime},t)\\}|r^{\text{sq.}},\mathfrak{B}\rangle\\\
&=(u_{k}^{\text{sq.}}(t)u_{k^{\prime}}^{\text{sq.}*}(t^{\prime})+u_{k^{\prime}}^{\text{sq.}}(t^{\prime})u_{k}^{\text{sq.}*}(t))\delta_{s,s^{\prime}}\delta_{\mathbf{k},-\mathbf{k}^{\prime}}\\\
\implies\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\\{\delta\hat{h}_{I}^{s}(\mathbf{k},t),\delta\hat{h}_{I}^{s^{\prime}}(\mathbf{k}^{\prime},t)\\}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}&=\delta_{s,s^{\prime}}\delta_{\mathbf{k}+\mathbf{k}^{\prime},0}\mathcal{Q}_{\delta
h}(t,t^{\prime},\mathbf{k})~{}.\end{split}$ (143)

It is important to note from the above equation that $\mathcal{Q}_{\delta
h}(t,t^{\prime},\mathbf{k})=(u_{k}^{\text{sq.}}(t)u_{k}^{\text{sq.}*}(t^{\prime})+u_{k}^{\text{sq.}}(t^{\prime})u_{k}^{\text{sq.}*}(t))$
which can be simplified as

$\begin{split}\mathcal{Q}_{\delta
h}(t,t^{\prime},\mathbf{k})&=2\Re(u_{k}^{\text{sq.}}(t)u_{k^{\prime}}^{\text{sq.}*}(t^{\prime}))=\frac{1}{k}\left(\cos\left(k(t-t^{\prime})\right)\cosh
2r_{k}-\cos\left(k(t+t^{\prime})-\phi_{k}\right)\sinh
2r_{k}\right)\end{split}$ (144)

which is eq.(90) from the main text of this paper. The non-squeezing case can
be reproduced from this result just by setting $r_{k}=0$ throughout the
analysis.

## References

  * [1] S. N. Bose, “Plancks Gesetz und Lichtquantenhypothese”, Zeitschrift für Physik 26 (1924) 178.