Quantum decoherence of phonons in Bose–Einstein condensates

The quantum field Hamiltonian for a rarefied, interacting, non-relativistic Bose gas is (see, e.g., [64])

$$\begin{eqnarray}\hat{H} & = & \displaystyle \int {\rm{d}}{\boldsymbol{r}}{\hat{\psi }}^{\dagger }({\boldsymbol{r}})\left[-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}^{2}+{ \mathcal V }({\boldsymbol{r}})\right]\hat{\psi }({\boldsymbol{r}})\\ & & +\displaystyle \frac{1}{2}\displaystyle \int {\rm{d}}{\boldsymbol{r}}{\rm{d}}{{\boldsymbol{r}}}^{{\prime} }{\hat{\psi }}^{\dagger }({\boldsymbol{r}}){\hat{\psi }}^{\dagger }({{\boldsymbol{r}}}^{{\prime} }){ \mathcal U }({{\boldsymbol{r}}}^{{\prime} }-{\boldsymbol{r}})\hat{\psi }({\boldsymbol{r}})\hat{\psi }({{\boldsymbol{r}}}^{{\prime} }),\end{eqnarray} \tag{ 1 }$$

where ${\hat{\psi }}^{\dagger }({\boldsymbol{r}})$ and $\hat{\psi }({\boldsymbol{r}})$ are the field operators creating and annihilating a bosonic atom at position ${\boldsymbol{r}};$ ${ \mathcal U }({\boldsymbol{r}})$ is the two-body potential; and ${ \mathcal V }({\boldsymbol{r}})$ is the trapping potential.

For clarity, we first consider a uniform gas occupying a volume V and later discuss how the analysis is straightforwardly extended to trapped systems. Although BECs can be created in a uniform trap [65], they are usually constrained in harmonic potentials. However, since studying uniform BECs is generally much simpler, these systems are often used as a first step in the theoretical analysis of trapped systems to describe certain properties through methods such as the local density approximation (see e.g., [66]) or by only probing small central sections of the trapped ultracold gas [67–69].

Considering a uniform system, we can substitute the plane-wave solutions $\hat{\psi }({\boldsymbol{r}})=\tfrac{1}{\sqrt{V}}{\sum }_{{\boldsymbol{k}}}{\hat{a}}_{{\boldsymbol{k}}}{{\rm{e}}}^{{\rm{i}}{\boldsymbol{k}}.{\boldsymbol{r}}}$ with ${\boldsymbol{k}}=2\pi {\boldsymbol{n}}/L;$ ${\boldsymbol{n}}=({n}_{x},{n}_{y},{n}_{z});$ and ${n}_{x},{n}_{y},{n}_{z}\in {\mathbb{N}}$ into (1) to obtain a complete description of the gas in terms of annihilation and creation operators in momentum space:

$$\begin{eqnarray}&&\hat{H}=\displaystyle \sum _{{\boldsymbol{k}}}\displaystyle \frac{{p}^{2}}{2m}{\hat{a}}_{{\boldsymbol{k}}}^{\dagger }{\hat{a}}_{{\boldsymbol{k}}}+\displaystyle \frac{1}{2V}\displaystyle \sum _{{\boldsymbol{k}},{{\boldsymbol{k}}}^{{\prime} },{\boldsymbol{q}}}{{ \mathcal U }}_{{\boldsymbol{q}}}{\hat{a}}_{{\boldsymbol{k}}+{\boldsymbol{q}}}^{\dagger }{\hat{a}}_{{{\boldsymbol{k}}}^{{\prime} }-{\boldsymbol{q}}}^{\dagger }{\hat{a}}_{{{\boldsymbol{k}}}^{{\prime} }}{\hat{a}}_{{\boldsymbol{k}}},\end{eqnarray} \tag{ 2 }$$

where ${{ \mathcal U }}_{{\boldsymbol{q}}}=\int { \mathcal U }({\boldsymbol{r}}){{\rm{e}}}^{-{\rm{i}}{\boldsymbol{q}}.{\boldsymbol{r}}}{\rm{d}}{\boldsymbol{r}}$.

Assuming the condensate to be macroscopically occupied, we next apply the Bogoliubov approximation [70]. This involves replacing ${\hat{a}}_{0}$ and ${\hat{a}}_{0}^{\dagger }$ with the c-number $\sqrt{N}$ and only retaining terms quadratic in ${a}_{{\boldsymbol{k}}\ne 0}$ and ${a}_{{\boldsymbol{k}}\ne 0}^{\dagger }$ (the higher-order terms are suppressed since they have fewer factors of $\sqrt{N}\gg 1$), where N is the number of atoms. In this approximation one also replaces the microscopic potential ${ \mathcal U }$ with an effective soft potential and expands the ${\boldsymbol{q}}=0$ component up to quadratic terms in the coupling constant $g=4\pi {{\hslash }}^{2}a/m$ where a is the s-wave scattering length and m is the atomic mass [64, 71].

The resulting Hamiltonian in ${a}_{{\boldsymbol{k}}\ne 0}$ and ${a}_{{\boldsymbol{k}}\ne 0}^{\dagger }$ can then be diagonalized by applying the following Bogoliubov transformation (see, e.g., [64])

$$\begin{eqnarray}&&{\hat{a}}_{{\boldsymbol{k}}}={u}_{k}{\hat{b}}_{{\boldsymbol{k}}}+{v}_{k}{\hat{b}}_{-{\boldsymbol{k}}}^{\dagger },\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\hat{a}}_{{\boldsymbol{k}}}^{\dagger }={u}_{k}{\hat{b}}_{{\boldsymbol{k}}}^{\dagger }+{v}_{k}{\hat{b}}_{-{\boldsymbol{k}}},\end{eqnarray} \tag{ 4 }$$

to obtain

$$\begin{eqnarray}&&\hat{H}={\epsilon }_{0}+\displaystyle \sum _{{\boldsymbol{k}}\ne 0}{\hslash }{\omega }_{k}{\hat{b}}_{{\boldsymbol{k}}}^{\dagger }{\hat{b}}_{{\boldsymbol{k}}},\end{eqnarray} \tag{ 5 }$$

where ${\hat{b}}_{{\boldsymbol{k}}}^{\dagger }$ and ${\hat{b}}_{{\boldsymbol{k}}}$ are the creation and annihilation operators for quasi-particles; ${\epsilon }_{0}$ is their ground state energy; c_s and ${\omega }_{k}$ are respectively the speed of sound and quasi-particle frequency

$$\begin{eqnarray}&&{({\hslash }{\omega }_{k})}^{2}:= {({c}_{s}{\hslash }k)}^{2}+{({{\hslash }}^{2}{k}^{2}/2m)}^{2},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{c}_{s}^{2}:= {gn}/m,\end{eqnarray} \tag{ 7 }$$

and u_k and v_k must satisfy

$$\begin{eqnarray}&&{u}_{k},{v}_{k}=\pm {\left(\displaystyle \frac{{{\hslash }}^{2}{k}^{2}/2m+{{mc}}_{s}^{2}\pm {\hslash }{\omega }_{k}}{2{\hslash }{\omega }_{k}}\right)}^{\tfrac{1}{2}},\end{eqnarray} \tag{ 8 }$$

where n is the number density of the gas.

As illustrated by (5), the Bose gas in this Bogoliubov approximation can be described by a non-interacting gas of quasi-particles, where the low momentum modes (${\hslash }k\ll {{mc}}_{s}$) are phonons travelling at speed c_s since ${\omega }_{k}\approx {c}_{s}k$ in this regime. Given that we have an ideal gas of quasi-particles, in thermal equilibrium the average occupation number ${N}_{{\boldsymbol{k}}}$ of quasi-particles carrying momentum ${\hslash }{\boldsymbol{k}}$ must satisfy

$$\begin{eqnarray}&&{N}_{{\boldsymbol{k}}}:= \langle {\hat{b}}_{{\boldsymbol{k}}}^{\dagger }{\hat{b}}_{{\boldsymbol{k}}}\rangle =\displaystyle \frac{1}{{{\rm{e}}}^{{\beta }_{k}}-1},\end{eqnarray} \tag{ 9 }$$

where ${\beta }_{k}:= {\hslash }{\omega }_{k}/{k}_{B}T$.

In the Bogoliubov approximation the phonons are non-interacting and so have infinite lifetimes. However, this approximation only keeps terms that are quadratic in ${\hat{a}}_{{\boldsymbol{k}}}$ and ${\hat{a}}_{{\boldsymbol{k}}}^{\dagger }$. If we instead also include the (more suppressed) cubic and quartic terms then, after the above Bogoliubov transformation (3), these terms will provide interactions between the quasi-particles, resulting in finite lifetimes.

Concentrating on just a single momentum mode ${\boldsymbol{q}}$ of the phonons, the cubic interaction terms for this mode are [47, 49, 56]

$$\begin{eqnarray}&&{\hat{H}}_{I}={\hat{b}}_{{\boldsymbol{q}}}{\hat{E}}_{{\boldsymbol{q}}}^{\dagger }+{\hat{b}}_{{\boldsymbol{q}}}^{\dagger }{\hat{E}}_{{\boldsymbol{q}}},\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}&&{\hat{E}}_{{\boldsymbol{q}}}:= \hat{A}+\hat{B}+\hat{L},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{\hat{A}}_{{\boldsymbol{q}}}^{\dagger }:= g\sqrt{\displaystyle \frac{n}{V}}\displaystyle \sum _{{\boldsymbol{k}},{{\boldsymbol{k}}}^{{\prime} }\ne \{{\bf{0}},{\boldsymbol{k}}\}}{{ \mathcal A }}_{k,{k}^{{\prime} }}{\hat{b}}_{{\boldsymbol{k}}}{\hat{b}}_{{{\boldsymbol{k}}}^{{\prime} }}{\delta }_{-{\boldsymbol{q}},{\boldsymbol{k}}+{{\boldsymbol{k}}}^{{\prime} }},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{\hat{B}}_{{\boldsymbol{q}}}^{\dagger }:= g\sqrt{\displaystyle \frac{n}{V}}\displaystyle \sum _{{\boldsymbol{k}},{{\boldsymbol{k}}}^{{\prime} }\ne \{{\bf{0}},{\boldsymbol{q}}\}}{{ \mathcal B }}_{k,{k}^{{\prime} }}{\hat{b}}_{{\boldsymbol{k}}}^{\dagger }{\hat{b}}_{{{\boldsymbol{k}}}^{{\prime} }}^{\dagger }{\delta }_{{\boldsymbol{q}},{\boldsymbol{k}}+{{\boldsymbol{k}}}^{{\prime} }},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\hat{L}}_{{\boldsymbol{q}}}^{\dagger }:= g\sqrt{\displaystyle \frac{n}{V}}\displaystyle \sum _{{\boldsymbol{k}},{{\boldsymbol{k}}}^{{\prime} }\ne \{{\bf{0}},{\boldsymbol{q}}\}}{{ \mathcal L }}_{k,{k}^{{\prime} }}{\hat{b}}_{{\boldsymbol{k}}}{\hat{b}}_{{{\boldsymbol{k}}}^{{\prime} }}^{\dagger }{\delta }_{{\boldsymbol{q}},{{\boldsymbol{k}}}^{{\prime} }-{\boldsymbol{k}}},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}{{ \mathcal A }}_{k,{k}^{{\prime} }} & := & \,{u}_{q}({v}_{k}{v}_{{k}^{{\prime} }}+{u}_{k}{v}_{{k}^{{\prime} }}+{v}_{k}{u}_{{k}^{{\prime} }})\\ & & +{v}_{q}({u}_{k}{v}_{{k}^{{\prime} }}+{v}_{k}{u}_{{k}^{{\prime} }}+{u}_{k}{u}_{{k}^{{\prime} }}),\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}{{ \mathcal B }}_{k,{k}^{{\prime} }} & := & \,{u}_{q}({u}_{k}{u}_{{k}^{{\prime} }}+{v}_{k}{u}_{{k}^{{\prime} }}+{u}_{k}{v}_{{k}^{{\prime} }})\\ & & +{v}_{q}({v}_{k}{v}_{{k}^{{\prime} }}+{v}_{k}{u}_{{k}^{{\prime} }}+{u}_{k}{v}_{{k}^{{\prime} }}),\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}{{ \mathcal L }}_{k,{k}^{{\prime} }}/2 & := & \,{u}_{q}({v}_{k}{u}_{{k}^{{\prime} }}+{u}_{k}{u}_{{k}^{{\prime} }}+{v}_{k}{v}_{{k}^{{\prime} }})\\ & & +{v}_{q}({u}_{k}{v}_{{k}^{{\prime} }}+{u}_{k}{u}_{{k}^{{\prime} }}+{v}_{k}{v}_{{k}^{{\prime} }}).\end{eqnarray} \tag{ 17 }$$

The resonant interactions ${\hat{b}}_{{\boldsymbol{q}}}{\hat{L}}_{{\boldsymbol{q}}}^{\dagger }$ and ${\hat{b}}_{{\boldsymbol{q}}}{\hat{B}}_{{\boldsymbol{q}}}^{\dagger }$ are the Landau and Beliaev interactions [40–55]. In the Landau process ${\hat{b}}_{{\boldsymbol{q}}}{\hat{L}}^{\dagger }$, a quasi-particle from mode ${\boldsymbol{q}}$ collides with a quasi-particle from another mode to create a higher-energy quasi-particle. Since this requires the thermal occupation of a quasi-particle mode, the process vanishes at zero temperature. On the other hand, in the Beliaev process ${\hat{b}}_{{\boldsymbol{q}}}{\hat{B}}^{\dagger }$, a quasi-particle of mode ${\boldsymbol{q}}$ spontaneously annihilates into two new quasi-particles with lower energies, which is analogous to parametric down-conversion in quantum optics [72] and can occur at absolute zero. Since the processes originate from (2), they can also be considered from the point-of-view of four-body interactions between the condensate and non-condensate atoms.

Treating the above single phonon mode ${\boldsymbol{q}}$ as an open quantum system, and the rest of the quasi-particle modes as its environment, we can decompose the Hamiltonian for the full system $\hat{H}$ in the following way:

$$\begin{eqnarray}&&\hat{H}={\hat{H}}_{S}+{\hat{H}}_{E}+{\hat{H}}_{I},\end{eqnarray} \tag{ 18 }$$

where ${\hat{H}}_{S}$ and ${\hat{H}}_{E}$ derive from (5), and respectively describe the free Hamiltonian of the considered single-mode phonon system and all other quasi-particle modes

$$\begin{eqnarray}&&{\hat{H}}_{S}:= {\hslash }{\omega }_{q}{\hat{b}}_{{\boldsymbol{q}}}^{\dagger }{\hat{b}}_{{\boldsymbol{q}}},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{\hat{H}}_{E}:= \displaystyle \sum _{{\boldsymbol{k}}\ne \{{\bf{0}},{\boldsymbol{k}}\}}{\hslash }{\omega }_{k}{\hat{b}}_{{\boldsymbol{k}}}^{\dagger }{\hat{b}}_{{\boldsymbol{k}}},\end{eqnarray} \tag{ 20 }$$

and ${\hat{H}}_{I}$ is defined in (10). We explicitly ignore the interaction terms between the states of the large system E, which describes all quasi-particle modes with momentum ${\hslash }{\boldsymbol{k}}\ne {\hslash }{\boldsymbol{q}}\ne 0$ and is assumed to be in thermal equilibrium.

To determine the evolution of the single-mode phonon system, we assume that the initial state of the full system is a product state $\hat{\rho }(0)={\hat{\rho }}_{S}(0)\otimes {\hat{\rho }}_{E}(0)$ and perform the Born–Markov approximation. That is, we assume that the coupling between E and S is weak, which is well justified for a rarefied Bose gas at low temperatures, and that the future evolution of ${\hat{\rho }}_{S}(t)$ does not depend on its past history. The latter is satisfied when the environment correlation time ${\tau }_{E}$ is much shorter than the time scale for significant change in S, which occurs when the environment is a large system maintained in thermal equilibrium, as expected here. The evolution of the single-mode system is then defined by the Lindblad master equation which, in diagonal form, is given by

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\hat{\rho }}_{S}}{{\rm{d}}{t}}=-\displaystyle \frac{{\rm{i}}}{{\hslash }}[{\hat{H}}_{S}^{{\prime} },{\hat{\rho }}_{S}]+\displaystyle \sum _{i=1}^{i=4}\left({\hat{c}}_{i}{\hat{\rho }}_{S}{\hat{c}}_{i}^{\dagger }-\displaystyle \frac{1}{2}\{{\hat{c}}_{i}^{\dagger }{\hat{c}}_{i},{\hat{\rho }}_{S}\}\right),\end{eqnarray} \tag{ 21 }$$

where ${\hat{H}}_{S}^{{\prime} }$ is the renormalized free Hamiltonian of the single-mode phonon system. Assuming the environment to be in thermal equilibrium, we have ${\hat{c}}_{1}=\sqrt{{\gamma }_{1}}{\hat{b}}_{{\boldsymbol{q}}}$, ${\hat{c}}_{2}=\sqrt{{\gamma }_{2}}{\hat{b}}_{{\boldsymbol{q}}}^{\dagger }$ where the rates ${\gamma }_{1}$ and ${\gamma }_{2}$ are related to environment correlation functions

$$\begin{eqnarray}&&{\gamma }_{1}=\displaystyle \frac{1}{{{\hslash }}^{2}}{\int }_{-\infty }^{\infty }{{\rm{d}}{t}}^{{\prime} }{{\rm{e}}}^{{\rm{i}}{\omega }_{q}{t}^{{\prime} }}\langle {\tilde{E}}_{{\boldsymbol{q}}}(t){\tilde{E}}_{{\boldsymbol{q}}}^{\dagger }(t-{t}^{{\prime} }){\rangle }_{E},\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{\gamma }_{2}=\displaystyle \frac{1}{{{\hslash }}^{2}}{\int }_{-\infty }^{\infty }{{\rm{d}}{t}}^{{\prime} }{{\rm{e}}}^{-{\rm{i}}{\omega }_{q}{t}^{{\prime} }}\langle {\tilde{E}}_{{\boldsymbol{q}}}^{\dagger }(t){\tilde{E}}_{{\boldsymbol{q}}}(t-{t}^{{\prime} }){\rangle }_{E},\end{eqnarray} \tag{ 23 }$$

and ${\tilde{E}}_{{\boldsymbol{q}}}$ is the operator defined in (11) but now in the interaction picture. The master equation (21) therefore describes the time evolution of a single-mode quasi-particle of a BEC in thermal equilibrium when taking into account Landau and Beliaev interactions (see also [56]).

By restricting our analysis to Gaussian states we can use the covariance matrix formalism to conveniently describe the dynamics of the single-mode phonon system. Such states are assumed in the relativistic quantum devices proposed in [26–28] that exploit the phonons of BECs, and can be readily created in BECs. For example, such states can be generated using Bragg spectroscopy [73–77] and in processes that are acoustic analogues to the dynamical Casimir effect (DCE) [78, 79].

The covariance matrix formalism is often used in quantum optics and continuous-variable quantum information science, and is a convenient mathematical framework in quantum phase space for describing Gaussian states and their dynamics. It utilizes the fact that Gaussian states are completely defined by just their first ${\boldsymbol{d}}$ and second ${\boldsymbol{\sigma }}$ statistical moments

$$\begin{eqnarray}&&{d}_{i}:= \langle {\hat{x}}_{i}\rangle \,=\,{\rm{Tr}}({\hat{x}}_{i}{\hat{\rho }}_{S}),\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}{\sigma }_{{ij}} & := & \displaystyle \frac{1}{2}\langle \{{\hat{x}}_{i},{\hat{x}}_{j}\}\rangle -\langle {\hat{x}}_{i}\rangle \langle {\hat{x}}_{j}\rangle \\ & \equiv & \displaystyle \frac{1}{2}{\rm{Tr}}(\{{\hat{x}}_{i}-{d}_{i},{\hat{x}}_{j}-{d}_{j}\}{\hat{\rho }}_{S}),\end{eqnarray} \tag{ 25 }$$

where ${\hat{x}}_{i}$ are the quadrature phase space operators which, for the single-mode system, are defined as

$$\begin{eqnarray*}&&\hat{{\boldsymbol{x}}}:= \left(\begin{array}{c}{\hat{x}}_{1}\\ {\hat{x}}_{2}\end{array}\right):= \displaystyle \frac{1}{2\kappa }\left(\begin{array}{cc}1 & 1\\ -{\rm{i}} & {\rm{i}}\end{array}\right)\left(\begin{array}{c}{\hat{b}}_{{\boldsymbol{q}}}\\ {\hat{b}}_{{\boldsymbol{q}}}^{\dagger }\end{array}\right).\end{eqnarray*}$$

These phase space operators then obey the commutation relations $[{\hat{x}}_{i},{\hat{x}}_{j}]=\tfrac{{\rm{i}}}{2{\kappa }^{2}}{{\rm{\Omega }}}_{{ij}}$ where ${\boldsymbol{\Omega }}$ is of symplectic form, and κ is a constant that is often taken to be $1/\sqrt{2}$ in quantum optics. By assuming Gaussian states we have, therefore, reduced the complete description of the state down from the infinite dimensional space of the density matrix to terms of the two-dimensional matrices ${\boldsymbol{\sigma }}$ and ${\boldsymbol{d}}$, which are referred to as the covariance matrix and displacement matrix of the state respectively.

To determine how the covariance and displacement matrices evolve for the single-mode phonon system, we first transform the Lindblad equation for the density matrix (21) into the phase space basis ${\hat{x}}_{i}$ and then use this equation in differential versions of (24) and (25) [80–82]. In fact it is possible to derive the time evolution equations for the covariance and displacement matrices of a general M-dimensional Gaussian state whose density matrix obeys (21): taking a general quadratic Hamiltonian⁵ ${\hat{H}}_{S}={H}_{0}+\kappa {\hat{{\boldsymbol{x}}}}^{T}{{\boldsymbol{H}}}_{1}+{\kappa }^{2}{\hat{{\boldsymbol{x}}}}^{T}{{\boldsymbol{H}}}_{2}\hat{{\boldsymbol{x}}}$ where H₀ is a constant; ${{\boldsymbol{H}}}_{1}$ is a $2M$-dimensional column-vector; and ${{\boldsymbol{H}}}_{2}$ is a $2M\times 2M$ real-symmetric matrix, it is straightforward to show that the covariance and displacement matrices evolve in the following way:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\boldsymbol{d}}}{{\rm{d}}{t}}={{\boldsymbol{ \mathcal H }}}_{1}+{\boldsymbol{A}}{\boldsymbol{d}},\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\boldsymbol{\sigma }}}{{\rm{d}}{t}}={\boldsymbol{A}}{\boldsymbol{\sigma }}+{\boldsymbol{\sigma }}{{\boldsymbol{A}}}^{{\boldsymbol{T}}}+{\boldsymbol{D}},\end{eqnarray} \tag{ 27 }$$

where ${{\boldsymbol{ \mathcal H }}}_{1}:= {\boldsymbol{\Omega }}{{\boldsymbol{H}}}_{1}/{\hslash }$. These equations have been derived previously in quantum optics for the case when ${{\boldsymbol{H}}}_{1}=0$ [80–82, 84]. In these quantum optics studies, the matrix ${\boldsymbol{A}}$ and symmetric matrix ${\boldsymbol{D}}$ are referred to as the drift and diffusion matrices respectively [82]. They are defined as

$$\begin{eqnarray}&&{\boldsymbol{D}}:= \displaystyle \frac{1}{4{\kappa }^{4}}{\boldsymbol{\Omega }}{\rm{Re}}({{\boldsymbol{C}}}^{\dagger }{\boldsymbol{C}}){{\boldsymbol{\Omega }}}^{{\boldsymbol{T}}},\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{\boldsymbol{A}}:= {{\boldsymbol{ \mathcal H }}}_{2}+{\boldsymbol{ \mathcal K }},\end{eqnarray} \tag{ 29 }$$

where

$$\begin{eqnarray}&&{\boldsymbol{ \mathcal K }}:= \displaystyle \frac{1}{2{\kappa }^{2}}{\boldsymbol{\Omega }}{\rm{Im}}({{\boldsymbol{C}}}^{{\boldsymbol{\dagger }}}{\boldsymbol{C}}),\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal H }}}_{2}:= \displaystyle \frac{1}{{\hslash }}{\boldsymbol{\Omega }}{{\boldsymbol{H}}}_{2},\end{eqnarray} \tag{ 31 }$$

and the matrix ${\boldsymbol{C}}$ is defined by ${\hat{c}}_{i}={C}_{{ij}}{\hat{x}}_{j}$.

The general solution of (26) is

$$\begin{eqnarray}&&{\boldsymbol{d}}(t)={\boldsymbol{X}}(t){{\boldsymbol{d}}}_{0}+{\boldsymbol{Y}}(t),\end{eqnarray} \tag{ 32 }$$

where, when ${\boldsymbol{A}}$ is independent of time,

$$\begin{eqnarray}&&{\boldsymbol{X}}(t)={{\rm{e}}}^{{\boldsymbol{A}}t},\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}{\boldsymbol{Y}}(t) & = & {\displaystyle \int }_{0}^{t}{{\rm{e}}}^{{\boldsymbol{A}}(t-s)}{\rm{d}}{s}{{\boldsymbol{H}}}_{1}\\ & = & {{\boldsymbol{A}}}^{-1}({{\rm{e}}}^{{\boldsymbol{A}}t}-1){{\boldsymbol{H}}}_{1},\end{eqnarray} \tag{ 34 }$$

so that the evolution of the displacement matrix is given by

$$\begin{eqnarray}&&{\boldsymbol{d}}(t)={{\rm{e}}}^{{\boldsymbol{A}}t}{{\boldsymbol{d}}}_{0}+{\int }_{0}^{t}{{\rm{e}}}^{{\boldsymbol{A}}(t-s)}{\rm{d}}{s}\,{\boldsymbol{H}}.\end{eqnarray} \tag{ 35 }$$

The equation of motion for the covariance matrix (27) is, in general, a time varying differential Lyapunov matrix equation, of which the general solution is

$$\begin{eqnarray}&&{\boldsymbol{\sigma }}(t)={\boldsymbol{X}}(t){{\boldsymbol{\sigma }}}_{0}{{\boldsymbol{X}}}^{T}(t)+{\boldsymbol{Z}}(t),\end{eqnarray} \tag{ 36 }$$

where

$$\begin{eqnarray}&&{\boldsymbol{X}}(t)={\boldsymbol{\Phi }}(t,0),\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&{\boldsymbol{Z}}(t)={\int }_{0}^{t}{\boldsymbol{\Phi }}(t,s){\boldsymbol{D}}{{\boldsymbol{\Phi }}}^{T}(t,s){\rm{d}}{s}.\end{eqnarray} \tag{ 38 }$$

When ${\boldsymbol{A}}$ is independent of time⁶ , ${\boldsymbol{\Phi }}(t,s)={{\rm{e}}}^{{\boldsymbol{A}}(t-s)}$ and so the solution to (27) can be written as⁷

$$\begin{eqnarray}&&{\boldsymbol{\sigma }}(t)={{\rm{e}}}^{{\boldsymbol{A}}t}{{\boldsymbol{\sigma }}}_{0}{{\rm{e}}}^{{{\boldsymbol{A}}}^{T}t}+{\int }_{0}^{t}{{\rm{e}}}^{{\boldsymbol{A}}(t-s)}{\boldsymbol{D}}{{\rm{e}}}^{{{\boldsymbol{A}}}^{T}(t-s)}{\rm{d}}{s},\end{eqnarray} \tag{ 39 }$$

for which an analytical expression can be obtained when ${\boldsymbol{A}}$ is diagonalizable [87].

A connection can be made with the evolution of the displacement vector and covariance matrix generated by a Gaussian unitary by neglecting all dissipative effects (${\boldsymbol{D}}=0$, ${\boldsymbol{A}}={{\boldsymbol{ \mathcal H }}}_{2}$). In this case, assuming for convenience that ${{\boldsymbol{H}}}_{2}$ and ${{\boldsymbol{H}}}_{1}$ are independent of time, the solutions (32) and (36) reduce to

$$\begin{eqnarray}&&{\boldsymbol{d}}(t)={\boldsymbol{S}}{{\boldsymbol{d}}}_{0}+{\boldsymbol{e}},\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&{\boldsymbol{\sigma }}(t)={\boldsymbol{S}}{{\boldsymbol{\sigma }}}_{0}{{\boldsymbol{S}}}^{T},\end{eqnarray} \tag{ 41 }$$

where ${\boldsymbol{S}}:= {{\rm{e}}}^{{{\boldsymbol{ \mathcal H }}}_{2}t}$ is the symplectic transformation corresponding to the free unitary evolution of the system, and ${\boldsymbol{e}}:= {\boldsymbol{B}}{{\boldsymbol{ \mathcal H }}}_{1}$ is a real vector with ${\boldsymbol{B}}:= {{\boldsymbol{ \mathcal H }}}_{2}^{-1}({{\rm{e}}}^{{{\boldsymbol{ \mathcal H }}}_{2}t}-1)$. The matrix ${{\boldsymbol{ \mathcal H }}}_{2}\,={\boldsymbol{\Omega }}{{\boldsymbol{H}}}_{2}/{\hslash }$ thus forms a symplectic algebra so that ${{\boldsymbol{ \mathcal H }}}_{2}$ is a (real) Hamiltonian matrix and ${{\boldsymbol{H}}}_{2}$ is a (real) symmetric matrix, which is also required by the Hermitian property of $\hat{H}$.

For our single-mode phonon system, ${\hat{H}}_{S}$ is given by (19) and thus ${H}_{0}=0$, ${{\boldsymbol{H}}}_{1}={\bf{0}}$ and ${{\boldsymbol{H}}}_{2}={\hslash }{\omega }_{q}^{{\prime} }{{\boldsymbol{I}}}_{2}$ where ${{\boldsymbol{I}}}_{2}$ is the two-dimensional identity matrix. The environment is also assumed to be in thermal equilibrium, and so the diffusion and drift matrices are given by

$$\begin{eqnarray}&&{\boldsymbol{D}}=\displaystyle \frac{1}{4{\kappa }^{2}}{\gamma }_{T}{{\boldsymbol{I}}}_{2},\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&{\boldsymbol{A}}=-\displaystyle \frac{1}{2}\gamma {{\boldsymbol{I}}}_{2}+{\omega }_{{\boldsymbol{q}}}^{{\prime} }{\boldsymbol{\Omega }},\end{eqnarray} \tag{ 43 }$$

where ${\gamma }_{T}:= {\gamma }_{1}+{\gamma }_{2};$ $\gamma := {\gamma }_{2}-{\gamma }_{1};$ and ${\omega }_{q}^{{\prime} }$ is the renormalized frequency of the single-mode phonon system.

Since the environment is in thermal equilibrium, the rates ${\gamma }_{1}$ and ${\gamma }_{2}$ are not independent but instead satisfy ${\gamma }_{1}={{\rm{e}}}^{{\beta }_{q}}{\gamma }_{2}$ [56, 88]. We can, therefore, write ${\gamma }_{T}$ as

$$\begin{eqnarray}&&{\gamma }_{T}:= \gamma \coth \left(\displaystyle \frac{1}{2}{\beta }_{q}\right)=\gamma (1+2{N}_{q}^{{\rm{th}}}),\end{eqnarray} \tag{ 44 }$$

where N^th_q is the average thermal occupation defined in (9). The matrix ${\boldsymbol{D}}$ is then given by

$$\begin{eqnarray}&&{\boldsymbol{D}}=\displaystyle \frac{(1+2{N}_{q})}{4{\kappa }^{2}}\gamma {{\boldsymbol{I}}}_{2}:= \gamma {{\boldsymbol{\sigma }}}_{\infty },\end{eqnarray} \tag{ 45 }$$

where ${{\boldsymbol{\sigma }}}_{\infty }$ is the covariance matrix of a single-mode thermal state.

Substituting the above drift and diffusion matrices for the single-mode phonon system into the general time-independent solutions (32) (using (33), (34)) and (39), the displacement vector and covariance matrix at time t are given by

$$\begin{eqnarray}&&{\boldsymbol{d}}(t)={{\rm{e}}}^{-\tfrac{1}{2}\gamma t}{\boldsymbol{R}}(t){{\boldsymbol{d}}}_{0},\end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray}&&{\boldsymbol{\sigma }}(t)={{\rm{e}}}^{-\gamma t}({\boldsymbol{R}}(t){{\boldsymbol{\sigma }}}_{0}{{\boldsymbol{R}}}^{T}(t))+(1-{{\rm{e}}}^{-\gamma t}){{\boldsymbol{\sigma }}}_{\infty },\end{eqnarray} \tag{ 47 }$$

where ${\boldsymbol{R}}(t):= {{\rm{e}}}^{\tfrac{1}{{\hslash }}{\boldsymbol{\Omega }}{{\boldsymbol{H}}}_{2}t}$ is the symplectic transformation corresponding to the free unitary evolution of the single-mode phonon system $U={{\rm{e}}}^{-\tfrac{{\rm{i}}}{{\hslash }}{H}_{S}t}$, and thus ${\boldsymbol{\Omega }}{{\boldsymbol{H}}}_{2}$ forms the symplectic algebra (it is a Hamiltonian matrix).

Since ${{\boldsymbol{H}}}_{2}={\hslash }{\omega }_{q}^{{\prime} }{{\boldsymbol{I}}}_{2}$ for our phonon system, ${\boldsymbol{R}}(t)\,=\cos ({\omega }_{q}t){{\boldsymbol{I}}}_{2}+\sin ({\omega }_{q}t){\boldsymbol{\Omega }}$, which is just the usual symplectic (and in this case rotational) transformation for the phase shift operator. From (47), due to the dissipative effects, this free evolution is now damped by ${{\rm{e}}}^{-\gamma t}$ and the state asymptotically approaches ${{\boldsymbol{\sigma }}}_{\infty }$:

$$\begin{eqnarray}&&{\boldsymbol{\sigma }}(t)={{\rm{e}}}^{-\gamma t}{{\boldsymbol{\sigma }}}_{0}+(1-{{\rm{e}}}^{-\gamma t}){{\boldsymbol{\sigma }}}_{\infty }.\end{eqnarray} \tag{ 48 }$$

This form of equation has been considered in quantum optics studies [89–91] but the rate γ in those studies is not the same as that for the phononic system studied here. This rate comes from the expressions for ${\gamma }_{1}$ and ${\gamma }_{2}$ (22)–(23) given in section 1.3. From these expressions, and assuming a continuum of modes, it can be easily shown that $\gamma ={\gamma }_{{\rm{B}}}+{\gamma }_{{\rm{L}}}$ where

$$\begin{eqnarray}&&{\gamma }_{{\rm{L}}}:= \displaystyle \frac{{g}^{2}n}{V{{\hslash }}^{2}}{\int }_{0}^{\infty }\pi {\rm{d}}{\omega }_{k}{p}_{k}{{ \mathcal L }}_{{kl}}^{2}\delta ({\omega }_{q}+{\omega }_{k}-{\omega }_{l})({N}_{l}^{{\rm{th}}}-{N}_{k}^{{\rm{th}}}),\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}&&{\gamma }_{{\rm{B}}}:= \displaystyle \frac{2{g}^{2}n}{V{{\hslash }}^{2}}{\displaystyle \int }_{0}^{\infty }\pi {\rm{d}}{\omega }_{k}{p}_{k}{{ \mathcal B }}_{{kl}}^{2}\delta ({\omega }_{q}-{\omega }_{k}-{\omega }_{l})\\ &&\quad \times \,(1+{N}_{l}^{{\rm{th}}}+{N}_{k}^{{\rm{th}}}),\end{eqnarray} \tag{ 50 }$$

which are just the Landau and Beliaev damping rates [40–56], where p_k is an assumed density of states. These rates have been explicitly calculated under various approximations. For example, in [49] the rates were calculated for a uniform three-dimensional BEC where it was found that, when ${k}_{B}T\ll {\hslash }{\omega }_{q}$ such that Beliaev damping dominates over Landau damping ${\gamma }_{{\rm{B}}}\gg {\gamma }_{{\rm{L}}}$, the total damping rate γ can be approximated as [49]

$$\begin{eqnarray}&&\gamma \approx {\gamma }_{{\rm{B}}}\approx \displaystyle \frac{3}{640\pi }\displaystyle \frac{{\hslash }{\omega }_{q}^{5}}{{{mnc}}_{s}^{5}}\left[1+{\left(\displaystyle \frac{{k}_{{B}}T}{{\hslash }{\omega }_{q}}\right)}^{3}\right],\end{eqnarray} \tag{ 51 }$$

which does not vanish at T = 0. On the other hand, in the opposite regime ${k}_{B}T\gg {\hslash }{\omega }_{q}$, Landau damping dominates over Beliaev damping, and for very high temperatures ${k}_{B}T\gg \mu \,\gg {\hslash }{\omega }_{q}$ (where $\mu ={gn}={{mc}}_{s}^{2}$ is the chemical potential), the total damping rate was found to be given by [49]

$$\begin{eqnarray}&&\gamma \approx {\gamma }_{{\rm{L}}}\approx \displaystyle \frac{3\pi }{8}\displaystyle \frac{{k}_{{B}}{Ta}}{{\hslash }{c}_{s}}{\omega }_{q},\end{eqnarray} \tag{ 52 }$$

whereas, for temperatures such that $\mu \gg {k}_{B}T\gg {\hslash }{\omega }_{q}$, the damping rate is found to be [49]

$$\begin{eqnarray}&&\gamma \approx {\gamma }_{{\rm{L}}}\approx \displaystyle \frac{3{\pi }^{3}}{40}\displaystyle \frac{{({k}_{B}T)}^{4}}{{mn}{{\hslash }}^{3}{c}_{s}^{5}}{\omega }_{q}.\end{eqnarray} \tag{ 53 }$$

Starting from the full Hamiltonian of a Bose gas (1), we have derived how the covariance matrix of a single-mode Gaussian phonon mode evolves in an isolated BEC (48). In the next section we investigate how this can then be used to find how quickly a prepared quantum state of phonons will decohere.

We now mention a few important steps in the above derivation of (48) from (1). The above discussion concentrated on a uniform BEC, such as that created in experiments performed in [65]. However, this can also be straightforwardly extended to more common trapped BECs, such as an harmonic trap. In this case we would start again from (1) but with ${ \mathcal V }({\boldsymbol{r}})$ replaced with the particular trapping potential, such as ${ \mathcal V }({\boldsymbol{r}})=\tfrac{1}{2}m({\omega }_{x}^{2}{x}^{2}+{\omega }_{y}^{2}{y}^{2}+{\omega }_{z}^{2}{z}^{2})$ for an anisotropic harmonic trap. The Bogoliubov transformations (3) can then be applied again to diagonalize the Hamiltonian in the Bogoliubov approximation (5), but now the coefficients u_k and v_k will be different to the uniform case (see e.g., [50, 51]). Going to next order in $\sqrt{N}$ we would then also find Beliaev and Landau interaction terms between the phonons as in the uniform case (10) but again with different coefficients [50, 51]. Therefore, for a non-uniform trapped system we will also end up with the evolution equation (48) for the covariance matrix of the phonons but with the damping rate of that particular trapped system. For example, the damping rate for non-uniform trapped systems has been analysed for low-energy excitations in [50, 51].

We also note that, even though the interaction Hamiltonian for our phonon system (10) is cubic in field modes, Gaussianity of the state will be preserved under the Born–Markov approximation. This is evident from the fact that the general solution to (27) has the form of the transformation brought about by a general Gaussian channel, which is a trace-preserving completely positive map that maps Gaussian trace-class operators onto Gaussian trace-class operators [92–102]. It is also illustrated by studying the Fokker–Planck equation for the Wigner function that derives from a general Markov master equation for the density operator, where it is found that the Wigner function continues to be Gaussian [89, 90]. The evolution of a Gaussian state in the Markov approximation is generally analysed for an interaction Hamiltonian that is bilinear in the field operators (see, e.g., [89, 90, 103]). However, the above analysis for a single-mode phonon state of a BEC clearly emphasizes that this is not a necessary condition for the preservation of Gaussianity. This is important since it means that, under realistic approximations, the simple description of the displacement and covariance matrices fully characterizing the state will persist throughout the state's evolution.

For a single-mode Gaussian state, the purity $\mu := {\rm{Tr}}({\rho }_{S}^{2})$ is given by

$$\begin{eqnarray}\mu & = & \displaystyle \frac{1}{4{\kappa }^{2}\sqrt{\det {\boldsymbol{\sigma }}}}\\ & = & \displaystyle \frac{1}{4{\kappa }^{2}s}\\ & = & \displaystyle \frac{1}{1+2{ \mathcal N }},\end{eqnarray} \tag{ 54 }$$

where ${ \mathcal N }$ is the so-called 'thermal' occupation of the state and s is its symplectic eigenvalue (the eigenvalue of the matrix $| {\rm{i}}{\boldsymbol{\Omega }}{\boldsymbol{\sigma }}|$). Using Williamson's theorem [107], any single-mode covariance matrix can be written in the general form [89, 90, 108]

$$\begin{eqnarray}&&{\boldsymbol{\sigma }}=\displaystyle \frac{1}{4{\kappa }^{2}\mu }\\ \quad \times \,\left(\begin{array}{cc}\cosh 2r+\sinh 2r\cos \psi & \sinh 2r\sin \psi \\ \sinh 2r\sin \psi & \cosh 2r-\sinh 2r\cos \psi \end{array}\right),\end{eqnarray} \tag{ 55 }$$

where r and ψ are defined by $\xi ={{r}{\rm{e}}}^{{\rm{i}}\psi }$, which is the squeezing parameter for a single-mode squeezing transformation. Using this general form of the covariance matrix in the solution (47) of its equation of motion, the purity of the single-mode state is found to evolve as [91, 105]

$$\begin{eqnarray}\mu (t) & = & {\mu }_{0}\left({{\rm{e}}}^{-2\gamma t}+\displaystyle \frac{{\mu }_{0}^{2}}{{\mu }_{\infty }^{2}}{(1-{{\rm{e}}}^{-\gamma t})}^{2}\right.\\ & & +{\left.\displaystyle \frac{2{\mu }_{0}}{{\mu }_{\infty }}{{\rm{e}}}^{-\gamma t}(1-{{\rm{e}}}^{-\gamma t})\cosh 2{r}_{0}\right)}^{-\tfrac{1}{2}},\end{eqnarray} \tag{ 56 }$$

where ${\mu }_{\infty }={(1+2{N}_{q}^{{\rm{th}}})}^{-1}=\tanh (\tfrac{1}{2}{\beta }_{q})$ is the purity of the thermal state ${{\boldsymbol{\sigma }}}_{\infty }$ that the covariance matrix asymptotically approaches. The purity will undergo a local minimum for squeezed states for which ${r}_{0}\gt \max [{\mu }_{0}/{\mu }_{\infty },{\mu }_{\infty }/{\mu }_{0}]$ at [91, 105, 106]

$$\begin{eqnarray}&&{t}_{\min }=\displaystyle \frac{1}{\gamma }\mathrm{ln}\left(\displaystyle \frac{\tfrac{{\mu }_{0}}{{\mu }_{\infty }}+\tfrac{{\mu }_{\infty }}{{\mu }_{0}}-2\cosh (2{r}_{0})}{\tfrac{{\mu }_{0}}{{\mu }_{\infty }}-\cosh (2{r}_{0})}\right),\end{eqnarray} \tag{ 57 }$$

which can provide a good characterization of the decoherence time of such states since this represents the point at which the state becomes most mixed, with any subsequent increase in the purity just reflecting the state being driven towards the state of the environment [91, 106]. To determine this decoherence time one just needs to know the initial purity and squeezing of the state; the temperature of the BEC (which defines ${\mu }_{\infty }$); and the damping rate γ. As an example, the quantum decoherence time for phonons in a particular BEC is calculated and presented in section 4. This BEC is that assumed in [27, 28] and the results can therefore be used to predict a quantum decoherence time for the phononic GW detector discussed in the Introduction.

An alternative characterization of the decoherence time can be provided by the nonclassical depth [109], which is a popular measure for quantifying the nonclassicality of a quantum state. This measure has the physical meaning of the number of thermal photons necessary to destroy the nonclassical nature of the quantum state [110]. For a general Gaussian state, the nonclassical depth τ detects the state as nonclassical if a canonical quadrature exists whose variance is below 1/2 [91] and, for a single-mode Gaussian state, it is given by

$$\begin{eqnarray}&&\tau =\max \left[\displaystyle \frac{1}{2}\left(1-\displaystyle \frac{{{\rm{e}}}^{-2r}}{\mu }\right),0\right].\end{eqnarray} \tag{ 58 }$$

From the time evolution of the covariance matrix (47), the nonclassical depth τ of the single-mode can be shown to evolve as [91]

$$\begin{eqnarray}&&\tau (t)=\displaystyle \frac{1}{2{\mu }_{\infty }}\left[{{\rm{e}}}^{-\gamma t}\left(1-\displaystyle \frac{{\mu }_{\infty }}{{\mu }_{0}}{{\rm{e}}}^{-2{r}_{0}}\right)+{\mu }_{\infty }-1\right].\end{eqnarray} \tag{ 59 }$$

Note that the time at which the state becomes 'classical' ($\tau =0$) is given by

$$\begin{eqnarray}&&{t}_{\tau =0}=\displaystyle \frac{1}{\gamma }\mathrm{ln}\left(\displaystyle \frac{1-\tfrac{{\mu }_{\infty }}{{\mu }_{0}}{{\rm{e}}}^{-2{r}_{0}}}{1-{\mu }_{\infty }}\right),\end{eqnarray} \tag{ 60 }$$

providing an alternative description for quantifying the decoherence time of quantum states.

As illustrated by the general covariance matrix of a single-mode Gaussian state (49), such states are fully defined by their first-statistical moment, purity, and squeezing $\xi ={{r}{\rm{e}}}^{{\rm{i}}\psi }$. For a single-mode state whose covariance matrix evolves as in (47), ψ is a constant of motion and r evolves as [91]

$$\begin{eqnarray}&&{\rm{\cosh }}(2r(t))=\mu (t)\left({{\rm{e}}}^{-\gamma t}\displaystyle \frac{{\rm{\cosh }}2{r}_{0}}{{\mu }_{0}}+\displaystyle \frac{1-{{\rm{e}}}^{-\gamma t}}{{\mu }_{\infty }}\right).\end{eqnarray} \tag{ 61 }$$

The average occupation of a single-mode state ${\boldsymbol{q}}$ is defined as ${N}_{{\boldsymbol{q}}}:= \langle {\hat{b}}_{{\boldsymbol{q}}}^{\dagger }{\hat{b}}_{{\boldsymbol{q}}}\rangle$. From the Lindblad master equation for the density operator of the single-mode state interacting with an environment in thermal equilibrium (21), it is easy to show that the average occupation of the single-mode phonon state evolves as

$$\begin{eqnarray}&&{N}_{{\boldsymbol{q}}}(t)={{\rm{e}}}^{-\gamma t}{N}_{{\boldsymbol{q}}}(0)+(1-{{\rm{e}}}^{-\gamma t}){N}_{q}^{{\rm{th}}},\end{eqnarray} \tag{ 62 }$$

where ${N}_{q}^{{\rm{th}}}:= {({{\rm{e}}}^{\beta {\omega }_{q}}-1)}^{-1}$ is the thermal occupation of the mode at temperature T. Since the average occupation is simply related to the average energy of the state, one can characterize the relaxation of the state from its evolution. Furthermore, since the average occupation is also defined by ${N}_{q}={\kappa }^{2}{\rm{Tr}}({\boldsymbol{\sigma }})+{\kappa }^{2}{{\boldsymbol{d}}}^{T}{\boldsymbol{d}}-1/2$, a single-mode Gaussian state at a time t is fully defined by its average occupation ${N}_{{\boldsymbol{q}}}$, purity μ and squeezing r.

It should be emphasized that the above estimate for the quantum decoherence the GW detector was calculated under various assumptions. In particular, a three-dimensional uniform BEC is assumed, whereas the GW detector was originally investigated in a one-dimensional setting. For one-dimensional Bose gases there is no condensation mechanism. Instead quasi one-dimensional BECs are investigated in experiments where the trapping frequencies in two of the dimensions are much greater than the third, making it possible to effectively integrate out the two dimensions to leave an approximate one-dimensional Gross–Pitaevskii equation for the condensate. Phonon damping has been investigated experimentally when moving from a three-dimensional setting to a quasi one-dimensional BEC [119, 120], where it was found that the damping rate was reduced, with the reduction depending on the temperature [120]. Attempts to theoretically model damping in quasi one-dimensional BECs are provided in [121–125].

A continuum of modes has also been assumed, which is an assumption that is often made in deriving phonon damping [40–56]. In particular, we have assumed that the phonon frequency of interest is ten times that of the fundamental frequency ${\omega }_{1}$ of the BEC, similar to that investigated in [27]. A discrete spectrum is likely to lead to an increase in decoherence times due to the reduction in phase space for the output states. For example, the discretization of the phonon spectrum in the radial directions of an elongated trap was studied in [125] to explain the experimental results of [120]¹⁰ . A large difference can be observed compared to continuous calculations for high temperatures, but the difference reduces at low temperatures.

The effect of the GW itself was also neglected in these calculations, which would be expected to extend the coherence life of the phonons of interest. How phonon interactions affect the squeezing created by the DCE, which is related to how the GWs create coherence of phonons [27], has been recently studied in [126].

Another important assumption has been the neglect of how three-body recombination of the BEC affects the phonon modes. For a three-dimensional BEC, three-body interactions cause the density of the BEC to evolve as [127]

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\rho (t)}{{\rm{d}}{t}}=-\gamma (t)\rho (t),\end{eqnarray} \tag{ 63 }$$

where $\gamma (t):= L{\rho }^{2}(t)$ and L has been found to be $5.8(1.9)\times {10}^{-42}\,{{\rm{m}}}^{6}\,{{\rm{s}}}^{-1}$ for a ${}^{87}\mathrm{Rb}$ BEC in the F = 1, ${m}_{f}=-1$ hyperfine state [128]. From (63), the half-life of the BEC is approximately $3/(2\gamma (0))$. Using ${c}_{s}=3.4\,\mathrm{mm}\,{{\rm{s}}}^{-1}$ we find that decay rate $\gamma (0)$ is less than the Beliaev damping rate ${\gamma }_{{\rm{B}}}$, and that the half-life of the BEC is over twice as long as the predicted time for full quantum decoherence of the phonons. Unlike for three-dimensional BECs, for which three-body recombination is constant at ultracold energies [129], in two and one-dimensional Bose gases three-body recombination can be vanishing at ultracold energies [130–132]. Therefore, moving to a quasi one- or two-dimensional BEC operating at low temperatures would be expected to reduce the three-body recombination rate [132].

As well as removing atoms from the trap, three-body recombination will also cause the BEC to heat up. In [133], the way in which three-body recombination, and in general n-body inelastic interactions, affects the phonons of a BEC was studied using an open system approach. A similar master equation is found to that of (21) but is derived using three-body inelastic interactions (and more generally n-body inelastic interactions), rather than the intrinsic two-body atomic interactions considered here, and the environment is the modes outside the trap. The relaxation rate was found to be comparable to the normal decay rate of three-body interactions, $\gamma (t)$ in (63), which is less than the Beliaev decay rate considered here for the three-dimensional detector at t = 0, as discussed above, and will continue to decrease as the particle number falls.

Of course, now that an estimate for the quantum decoherence time has been calculated for this frequency range we can also begin to investigate how the detector could be modified in order to increase this time. For example, one option might be to increase ${\omega }_{1}$ such that the initial state involves a lower energy mode of the cavity for which Beliaev damping will be suppressed. However, this would require modifications to c_s and the trap geometry, such as increasing the speed of sound and reducing the effective size of the trap, which would lead to more rapid three-body recombination. Other options to increase the quantum decoherence time could be to squeeze the environment [91], use larger mass BECs such as ${}^{174}\mathrm{Yb}$, and to utilise lower dimensional traps as described above.