%~Mouliné par MaN_auto v.0.32.4 2024-11-12 12:11:29
\documentclass[CRPHYS, Unicode, XML, screen, biblatex]{cedram}

\addbibresource{crphys20240218.bib}


\usepackage{xcolor}
\usepackage{soul}



\newcommand{\br}{\mathbf{r}}
\newcommand{\clE}{\mathcal{E}}
\newcommand{\clL}{\mathcal{L}}
\newcommand{\rconst}{\mathrm{const}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

%\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\let\oldhat\hat
\renewcommand*{\hat}[1]{\mathchoice{\widehat{#1}}{\widehat{#1}}{\oldhat{#1}}{\oldhat{#1}}}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\makeatletter\def\TITREspecial{relax}\def\cdr@specialtitle@english{Simulating gravitational problems with condensed matter analog models: a special issue in memory of Renaud Parentani (1962-2020)}
\def\cdr@specialtitle@french{Simuler des problèmes gravitationnels avec des modèles analogues en matière condensée : un numéro spécial en mémoire de Renaud Parentani (1962-2020)}
\makeatother{}


\title{The Dynamical Casimir Effect in quasi-one-dimensional Bose condensates: the breathing ring}

\alttitle{L'effet Casimir dynamique dans les condensats de Bose quasi-unidimensionnels : l'anneau de respiration}

\author{\firstname{Manuele} \lastname{Tettamanti}\CDRorcid{0000-0002-0889-9060}}
\address{Dipartimento di Fisica ``Giuseppe Occhialini'', Universit\`{a} di Milano-Bicocca and INFN - Sezione di Milano-Bicocca - Piazza della Scienza 3, 20126 Milano, Italy}


\author{\firstname{Alberto} \lastname{Parola}\CDRorcid{0000-0001-7844-6916}\IsCorresp}
\address{Dipartimento di Scienza e Alta Tecnologia and To.Sca.Lab, Universit\`{a} degli Studi dell'Insubria - Via Valleggio 11, 22100 Como, Italy}
\email[A. Parola]{alberto.parola@uninsubria.it}

\keywords{Casimir Effect, Bose Einstein Condensates, vacuum fluctuations}

\altkeywords{Effet Casimir, Condensats Bose Einstein, fluctuations du vide}

\begin{abstract}
We present a detailed investigation of one of the cleanest examples where it is possible to detect the ``analog'' Dynamical Casimir Effect in a Bose--Einstein condensate: an ultracold atom gas in toroidal confinement. The analytical solution of the time dependent Gross--Pitaevskii equation allows to follow the time evolution of the phonon spectrum and shows that periodic oscillations of the ring radius do not induce modulations in the density profile but give rise to the mixing of clockwise and anticlockwise modes, leading to the creation of pairs of entangled phonons in a squeezed vacuum state, if the drive frequency equals twice the frequency of the phonon mode. The Dynamical Casimir Effect is predicted to occur in the weakly interacting regime, where the Gross--Pitaevskii equation provides a faithful description of the many body dynamics. In the strong coupling limit, when the ultracold gas behaves as hard core bosons, the effect disappears and no amplification occurs. The presence of symmetry-breaking perturbations and finite temperature effects are also considered, as well as the comparison with the classical phenomenon of parametric amplification.
\end{abstract}

\begin{altabstract}
Nous présentons une enquête détaillée sur l'un des exemples les plus clairs où il est possible de détecter l'effet Casimir Dynamique \og analogue\fg{} dans un condensat de Bose--Einstein : un gaz d'atomes ultrafroids en confinement toro\"idal. La solution analytique de l'équation de Gross--Pitaevskii dépendant du temps permet de suivre l'évolution temporelle du spectre des phonons et montre que des oscillations périodiques du rayon de l'anneau n'induisent pas de modulations dans le profil de densité mais donnent lieu au mélange des modes dans le sens horaire et antihoraire, conduisant à la création de paires de phonons intriqués dans un état de vide comprimé, si la fréquence de l'entraînement est égale à deux fois la fréquence du mode de phonon. L'effet Casimir Dynamique est prévu se produire dans le régime d'interaction faible, où l'équation de Gross--Pitaevskii fournit une description fidèle de la dynamique à plusieurs corps. Dans la limite du couplage fort, lorsque le gaz ultrafroid se comporte comme des bosons de cœur dur, l'effet disparaît et aucune amplification ne se produit. La présence de perturbations de rupture de symétrie et des effets de température finie sont également considérés, ainsi que la comparaison avec le phénomène classique d'amplification paramétrique.
\end{altabstract}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
\maketitle

\section{Introduction}

The non-trivial nature of the vacuum marks a crucial difference between the classical and the quantum world~\cite{milonni}. Unfortunately, a direct evidence of the quantum fluctuations of the vacuum state is difficult to achieve unless some mechanism leads to their amplification. The Dynamical Casimir Effect~\cite{moore} might prove instrumental to reach this goal and considerable efforts were devoted to its detection in superconducting devices~\cite{lath,wilson}, while other intriguing effects related to vacuum fluctuations, like the Hawking~\cite{hawking} and Unruh~\cite{unruh} effects, have not been observed so far~\cite{nori}. In order to circumvent the difficulties in the direct experimental measurement of vacuum fluctuations, the idea of transferring the same mechanism into a different context has been widely explored. In this respect, ultracold atoms represent a privileged testing ground because they allow to achieve extremely low temperatures and can be efficiently manipulated by optical techniques. In fact, the detection of ``analog'' Dynamical Casimir Effect~\cite{jaskula} and the evidence of ``analog'' Hawking radiation~\cite{stein16,stein18} in ultracold atoms have been recently claimed, although the interpretation of the results always requires considerable care~\cite{carusotto,jacob}, especially when these subtle quantum mechanical effects are embedded in a condensed matter environment~\cite{epl,prd}.

A recent proposal~\cite{demler} suggested to investigate the ``analog'' Dynamical Casimir Effect (DCE) in ultracold atom gases by confining a Bose--Einstein Condensate (BEC) in a flat quasi-one dimensional trap (box) whose walls undergo periodic oscillations. The setting is strongly reminiscent of the original electromagnetic set-up. However, in this case the condition that the wall peak velocity should be comparable to the velocity of light is replaced by the analog condition on the much smaller sound velocity, which makes the experiment feasible. The system is expected to display the known phenomenon of parametric resonance, leading to an amplification of the virtual phonons stemming from vacuum fluctuations. While parametric resonance is ubiquitous in physics, from swings to electric circuits, the spontaneous DCE is much more subtle, being an intrinsically quantum phenomenon involving the amplification of \emph{vacuum fluctuations} and the creation of entangled pairs of the field quanta. In the case of the electromagnetic field, DCE is realized through the periodic motion of the mirrors confining an empty cavity in the absence of radiation: Virtual photons become real and the initial trivial vacuum state evolves into a squeezed vacuum state. It is crucial that the moving mirrors are uncharged and therefore do not directly couple to the electromagnetic field, otherwise the accelerated charges would emit photons, hiding the DCE. In the ``analog'' case of a BEC the corresponding requirement is that the phonon field must be in the vacuum state while the BEC is externally manipulated. Such a condition is very hard to achieve, because as soon as we mechanically interact with a classical or quantum fluid, sound waves are unavoidably generated. Indeed, the first experiment addressing this problem detected correlated phonon pairs traveling in opposite directions after squeezing the condensate~\cite{jaskula}, while in the alternative set-up~\cite{demler} the BEC is confined in a flat trap whose walls are periodically shaken. In both cases, the mechanical excitation of the BEC gives rise to waves in the fluid, which do not rely upon the quantum nature of the state. Therefore it is not clear whether the emerging phonons are going to display the highly non trivial entanglement features characterizing the DCE. A further, fundamental difference between the DCE in the electromagnetic framework and the analog effect in BEC's is that in the former case the field quanta (photons) are directly accessible in experiments while in the latter, phonons are just the elementary excitations of the condensate, i.e. quasiparticles, whose attributes may be inferred only by measuring the properties of the atoms in the BEC.

In this paper we investigate an alternative geometrical configuration where these problem are, in principle, overcome. The BEC is confined in a narrow toroidal trap, whose radius periodically oscillates with an external frequency $\Omega$ and very small amplitude: Symmetry is not broken and then the condensate remains uniform during the oscillations. If $\Omega$ equals twice a quantized eigenfrequency of the BEC, phonon pairs of opposite momenta are created on a secular timescale proportional to the inverse of the (small) amplitude of the oscillations. We show that if initially the BEC is in its ground state, the evolution indeed gives rise to a squeezed vacuum state characterized by an exponentially growing population of the two resonant modes. As already noticed, a hallmark of this squeezed state is the presence of strong correlations in the velocity-velocity correlation function of the atoms in the condensate~\cite{iac}. A unique feature of this two-mode squeezed state is the absence of fluctuations in the difference between the population of the two resonant modes: While the condensate displays quantum fluctuations in the number of atoms on each of the two modes, in an ideal experiment, the difference between the populations of the two modes should have vanishing fluctuations.

The paper is organized as follows: in Section~\ref{sect2} we first describe the set-up we are investigating and the reduction to a one dimensional effective model. We discuss both the weak coupling case, accurately reproduced by the cubic Gross--Pitaevskii equation, and the strong coupling, Tonks--Girardeau, limit. In Section~\ref{sect3} the properties of the ground state and of the excitation spectrum are reviewed. In Section~\ref{sect4} a few experimentally accessible properties of the condensate are examined. In Section~\ref{sect-time} we investigate the effects of a small periodic oscillation of the radius of the toroidal trap providing an explicit solution of the Gross--Pitaevskii equation in the resonant case. Section~\ref{sect6} illustrates the modifications in the velocity correlations induced by resonant oscillations, offering an interpretation in terms of analog DCE. The results of our quantum model are compared with a purely classical treatment in Section~\ref{sect7}, while Section~\ref{sect8} contains a brief summary of the main results and some final considerations.

\section{BEC in a breathing ring}\label{sect2}
We investigate a BEC in a narrow toroidal trap, confined by a time dependent harmonic potential of the form:
\begin{equation}\label{trap}
V_t(r,z) = \frac{1}{2} M \left [\omega_r^2 \,(r-R_t)^2 + \omega^2_z \,z^2\right]
\end{equation}
where $(r,\theta,z)$ are cylindrical coordinates and $R_t$ is the average radius which may depend on time. $M$ is the mass of each boson and $(\omega_r,\omega_z)$ the trap frequencies along $(r,z)$ respectively. Traps of toroidal geometry were already studied and realized in an atomic condensates~\cite{tori3,tori2,toro,tori4,tori1}. Although the details of the trapping potential depend on the experimental set-up, in the limit of strong confinement, the overall shape can be always parametrized by a form similar to Eq.~\eqref{trap}.

At weak coupling, the dynamics of the BEC is governed by the usual three dimensional Gross--Pitaevskii equation (GPE)\cite{stringa}:
\begin{equation}
i\hbar \frac{\partial\psi_t(\br)}{\partial t} = -\frac{\hbar^2}{2M} \nabla^2 \psi(\br) + V_t(\br)\,\psi_t(\br) + g_{3D} \,|\psi_t(\br)|^2 \psi_t(\br)
\end{equation}
where $g_{3D} =\frac{4\pi\hbar^2 a_s}{M} $ is expressed in terms of the (positive) scattering length $a_s$. The condensate wavefunction is normalized to the total number of particles: $\int d\br \,|\psi_t(\br)|^2 = N_0$. In the case of a narrow trap it is customary to simplify the problem by dimensional reduction~\cite{sala}: the solution is approximately represented as the product of the solution $G_t(r,z)$ of the transverse Schr\"odinger equation:
\begin{equation}
i\hbar\frac{\partial G_t(r,z)}{\partial t} = -\frac{\hbar^2}{2M} \nabla_\perp^2 G_t(r,z) + V_t(r,z)\,G_t(r,z)
\end{equation}
times an unknown periodic function $\phi_t(\theta)$ of the polar angle $\theta$: $\psi_t(\br) = \phi_t(\theta)\,G_t(r,z)$. The unknown angular wavefunction $\phi_t(\theta)$ is correctly normalized $\int_0^{2\pi} d\theta \,|\phi_t(\theta)|^2 = N_0$ provided $\int_0^\infty dr\,r\int_{-\infty}^\infty dz \,|G_t(r,z)|^2 = 1$. By substituting this parametrization in the GPE, multiplying by $G_t^*(r,z)$ and integrating over the $(r,z)$ coordinates, we recover the usual one dimensional effective GPE \footnote{The ``Hubble friction term'', found in Ref.~\cite{jacobson} is not present in our formalism due to the different choice of variables, see also Section~7.}:
\begin{equation}\label{gpe1d}
i\hbar \frac{\partial\phi_t(\theta)}{\partial t} = -\frac{\hbar^2}{2MR_t^2} \frac{\partial^2\phi_t(\theta)}{\partial\theta^2} +\frac{g}{R_t} \,|\phi_t(\theta)|^2 \phi_t(\theta)
\end{equation}
where $g = \frac{g_{3D}}{2\pi\sigma_r\sigma_z}$ with $\sigma_{r,z} = \sqrt{\frac{\hbar}{M\omega_{r,z}}}$. In order to obtain this expression we have approximated
\[
\int_0^\infty dr\,r\int_{-\infty}^\infty dz \,\frac{|G_t(r,z)|^2}{r^2} \sim R_t^{-2}\quad \text{and}\quad \int_0^\infty dr\,r\int_{-\infty}^\infty dz \,|G_t(r,z)|^4 \sim \left[2\pi\sigma_z\sigma_r\,R_t\right]^{-1}.
\]
These estimates are asymptotically exact for narrow traps $R_t \gg \sigma_r$, showing that the time dependence of the trap radius $R_t$ affects both the kinetic and the interaction terms. Note that this effective GPE equation is compatible with the periodicity of $\phi_t(\theta)$ in the interval $\theta\in[0,2\pi]$. Remarkably, if instead of the usual cubic GPE we start from the quintic GPE~\cite{menotti,sala1d} appropriate at strong coupling, the same argument shows that the effect of a time dependent radius is just to rescale the hamiltonian by a factor $(R_0/R_t)^2$ because both the kinetic and the interaction term scale in the same way, implying that the time evolution becomes trivial in this limit. The same conclusion can be reached by considering the Tonks--Girardeau (TG) fluid~\cite{girardeau}, i.e. a fluid of hard core bosons in one dimension, where the condensate is mapped into a free Fermi gas. Each single particle orbital $\phi_t(\theta)$ evolves in time according to the one dimensional Schr\"odinger equation:
\begin{equation}
i\hbar \frac{\partial\phi_t(\theta)}{\partial t} = -\frac{\hbar^2}{2MR_t^2} \frac{\partial^2\phi_t(\theta)}{\partial\theta^2}
\end{equation}
If at $t=0$ the system is set in the ground state, characterized by a set of occupied free particle states $\frac{e^{i m\theta}}{\sqrt{2\pi}}$, the time evolution just changes the amplitude of each single particle state by a time dependent phase factor. Therefore, the physical properties of the state do not change under the external force: no resonance is present. So we expect that, in this set-up, the Dynamical Casimir Effect is absent in the strongly interacting, hard-core limit, in agreement with the conclusion reached by the analysis of the \emph{quintic} GPE, as opposed to the usual cubic one. More generally, following Ref.~\cite{sala1d}, we can write a non-linear Schr\"odinger equation able to describe the dynamics of a condensate at intermediate couplings interpolating between the weakly interacting limit (GPE) and the strong coupling regime (TG).

Finally we note that an effective one dimensional equation of the same form as~\eqref{gpe1d} is recovered by fixing $R_t=R_0$ and varying the trapping frequencies, e.g. $\omega_z(t)$~\cite{butera}. In this case the kinetic term is time independent while the effective one dimensional coupling constant changes in time: $g(t)\propto \sqrt{\omega_z(t)}$. Also in this set-up the cylindrical symmetry of the problem is preserved and also in this case the effects of the perturbation disappear at strong coupling, when $g(t)\to\infty$.

\section{Ground state and excitation spectrum}\label{sect3}

The one dimensional GPE can be formally derived by use of the second quantization formalism\cite{fetter}, starting from the Grand Canonical Hamiltonian operator describing a gas of bosonic particles on a ring with contact interactions:
\begin{equation}\label{hamw}
\hat H = \int_0^{2\pi} d\theta \left\{\hat\psi^\dagger(\theta) \left [-\frac{\hbar^2}{2MR_t^2} \frac{d^2}{d\theta^2} -\mu\right] \hat\psi(\theta) + \frac{g}{2R_t}\left [\hat\psi^\dagger(\theta)\right]^2\left [\hat\psi(\theta)\right]^2 \right\}
\end{equation}
where $\hat \psi(\theta)$ are field operators obeying the canonical commutation relations:
\begin{equation}\label{commuta}
\left [\hat\psi(\theta),\hat\psi^\dagger(\theta^\prime)\right] = \delta(\theta-\theta^\prime)
\end{equation}
and $\mu$ is the chemical potential.

Let us start by considering a time independent trapping potential: $R_t=R_0$. In order to put~\eqref{hamw} in diagonal form we write the field operator as
\begin{equation}\label{split0}
\hat\psi(\theta) = \Phi(\theta) +\delta\hat\psi(\theta)
\end{equation}
where the ``condensate wavefunction'' $\Phi(\theta)$ is a suitable complex function. Clearly, the canonical commutation relations are satisfied for any choice of $\Phi(\theta)$ provided $\delta\hat\psi(\theta)$ obeys Eq.~\eqref{commuta}. Substituting into~\eqref{hamw} and neglecting terms of higher order in the fluctuations, we get
\begin{equation}\label{sum}
\hat H = E_0 + \hat H_1 + \hat H_2
\end{equation}
where
\begin{align}
E_0 &=
\int_0^{2\pi} d\theta \, \Phi^*(\theta) \Big (-\frac{\hbar^2}{2MR_0^2} \frac{d^2}{d\theta^2}-\mu + \frac{g}{2R_0} \,|\Phi(\theta)|^2 \Big) \Phi(\theta) \\
\hat H_1 &=
\int_0^{2\pi} d\theta \,\,\delta\hat\psi^\dagger(\theta)\,
\Big (-\frac{\hbar^2}{2MR_0^2} \frac{d^2}{d\theta^2} -\mu + \frac{g}{R_0}\,|\Phi(\theta)|^2 \Big)\Phi(\theta) + h.c.
\end{align}
\begin{multline}\label{h2}
\hat H_2 = \int_0^{2\pi} d\theta \Bigg \{\delta\hat\psi^\dagger(\theta) \Big (-\frac{\hbar^2}{2MR_0^2} \frac{d^2}{d\theta^2} + 2\frac{g}{R_0}\,|\Phi(\theta)|^2 -\mu\Big)\delta\hat\psi(\theta)\\
 +\, \frac{g}{2R_0} \, \left (\Phi(\theta)^2 \, \delta\hat\psi^\dagger(\theta)^2 + \Phi^*(\theta)^2 \,
\delta\hat\psi(\theta)^2 \right) \Bigg \}\quad
\end{multline}
Setting $\hat H_1=0$ we obtain the stationary Gross--Pitaevskii equation for the condensate wavefunction $\Phi(\theta)$. The chemical potential $\mu$ is defined by the normalization condition $\int_0^{2\pi} d\theta \,|\Phi(\theta)|^2 = N_0$. In our case, the solution is trivial because of rotational invariance: 
\[
\Phi(\theta) = \sqrt{\frac{N_0}{2\pi}}\quad\text{and}\quad\mu=\frac{N_0\,g}{2\pi R_0}=g\,\rho_0,
\]
where $\rho_0$ is the one dimensional number density of bosons.

Having determined the condensate wavefunction, let us investigate the shape of the excitation spectrum by solving the eigenvalue problem for $\hat H_2$. Following Bogoliubov~\cite{stringa}, we express the field operator in terms of creation and annihilation operators by
\begin{equation}\label{expa}
\delta\hat\psi(\theta) = \sum_n \left \{ D_n(\theta)\,\hat b_n + E_n^*(\theta)\hat b^\dagger_n\right \}
\end{equation}
where $D_n(\theta)$ are complex functions to be determined, $n$ labels the eigenmodes of the system and $\hat b_n$ are bosonic annihilation operators satisfying $[\hat b_n,\hat b^\dagger_m] = \delta_{n,m}$. The canonical commutation relation~\eqref{commuta} are correctly fulfilled provided the ortho-normalization conditions hold:
\begin{align}\label{norma}
\begin{split}
\sum_n \left \{ D_n(\theta) E^*_n(\theta^\prime) - E_n^*(\theta) D_n(\theta^\prime)
\right \} &= 0\\
\sum_n \left \{ D_n(\theta) D^*_n(\theta^\prime) - E_n^*(\theta) E_n(\theta^\prime)
\right \} &= \delta(\theta-\theta^\prime)
\end{split}
\end{align}
Inverting, we can write
\begin{equation}\label{inverse}
\hat b_n = \int_0^{2\pi} d\theta\, \left \{ D^*_n(\theta)\,\delta\hat\psi(\theta) - E^*_n(\theta) \delta\hat\psi(\theta)^\dagger\right \}
\end{equation}
while the commutation relations imply:
\begin{align}
\int_0^{2\pi} d\theta\, \left \{ D_n(\theta) E_{n^\prime}(\theta) - E_n(\theta) D_{n^\prime}(\theta)
\right \} &= 0 \label{norma1} \\
\int_0^{2\pi} d\theta\, \left \{ D^*_n(\theta) D_{n^\prime}(\theta) - E^*_n(\theta) E_{n^\prime}(\theta)
\right \} &= \delta_{n,n^\prime}\label{norma2}
\end{align}
The aim of the Bogoliubov transformation is to express the quadratic Hamiltonian $\hat H_2$~\eqref{h2} in diagonal form
\begin{equation}\label{diag}
\hat H_2 = \sum_n \,\hbar \omega_n \,\hat b^\dagger_n \hat b_n + \clE_0
\end{equation}
where $\omega_n$ are real quantities. The easiest way to determine the coefficients $D_n(\theta),E_n(\theta)$ and the excitation frequencies $\omega_n$ is to take the commutator of $\hat H_2$ with the field operator $\delta\hat \psi(\theta)$. The resulting eigenvalue system is
\begin{equation}\label{bogo}
\begin{aligned}
h_2^0 D_n(\theta) + g\,\rho_0 E_n(\theta) &= \hbar\omega_n \,D_n(\theta) \\
h_2^0 E_n(\theta) + g\,\rho_0 D_n(\theta) &= -\hbar\omega_n \,E_n(\theta)
\end{aligned}
\end{equation}
with
\begin{equation}\label{h20}
h_2^0= -\frac{\hbar^2}{2MR_0^2} \frac{d^2}{d\theta^2} +g\,\rho_0
\end{equation}
This eigensystem defines the set of solutions $D_n(\theta),E_n(\theta)$ and the corresponding eigenfrequencies $\omega_n$ apart from the normalization condition~\eqref{norma2} which must be imposed a posteriori. Note that a stable equilibrium solution must allow only for positive eigenfrequencies. It is clear from Eqs.~\eqref{bogo} that if a solution $(\omega,D,E)$ exists, formally, another ``negative norm'' (unphysical in this framework) solution $(-\omega,E^*,D^*)$ is also present. Note also the presence of a vanishing norm, trivial eigenvalue with $\omega=0$ with eigenfunctions $D(\theta)=-E(\theta)=\rconst$.

The vacuum fluctuation energy $\clE_0$ can be easily evaluated by taking the expectation value of~\eqref{diag} on the vacuum of the $\delta\hat\psi(\theta)$ operators. This gives:
\begin{equation}
\clE_0 = - \sum_n\,\hbar\omega_n\,\int_0^{2\pi} d\theta\, |E_n(\theta)|^2
\end{equation}

The solution of the Bogoliubov equations~\eqref{bogo} in a uniform system is:
\begin{equation}\label{bogo2}
D_n(\theta) = D_n \,\frac{e^{in\theta}}{\sqrt{2\pi}}\qquad\qquad E_n(\theta) = E_n \frac{\,e^{in\theta}}{\sqrt{2\pi}}
\end{equation}
where $n=0,\pm 1,\pm 2\dots$ and the amplitudes $(D_n,E_n)$ satisfy the algebraic eigenvalue system:
\begin{equation}\label{eig}
\begin{bmatrix}
\frac{\hbar^2 n^2}{2MR_0^2} + g\rho_0 & g\,\rho_0 \\
-g\,\rho_0 & -\frac{\hbar^2 n^2}{2MR_0^2} - g\rho_0
\end{bmatrix}
\begin{bmatrix}
D_n^\eta\\
E_n^\eta
\end{bmatrix} = \eta\,\hbar\omega_n\,
\begin{bmatrix}
D_n^\eta\\
E_n^\eta
\end{bmatrix}
\end{equation}
where the label $\eta=\pm$ identifies the two branches of positive and negative frequency. The secular equation provides the usual dispersion relation
\begin{equation}\label{omega}
\hbar\omega_n = \sqrt{\frac{g\rho_0\hbar^2 n^2}{MR_0^2} + \left [\frac{\hbar^2 n^2}{2MR_0^2}\right]^2}
\end{equation}
All the eigenvalues are real implying dynamical stability of the GPE equation~\cite{stab}. As stated in Eq.~\eqref{eig}, the secular equation actually has two solutions for each $n\ne 0$: $\pm \hbar\omega_n$. If $n=0$ the only solution has vanishing eigenvalue and vanishing norm and therefore does not represent a physical eigenmode of the system. The physical branch of the eigenvalues corresponds only to the positive sign $\eta=+$, because the negative frequency solutions give rise to eigenvectors of negative norm (thereby not satisfying the wanted commutation relations). In the following, having discarded the negative frequency branch, we drop the label $+$ on the eigenvectors. Finally, the symmetry $\theta \to -\theta$ implies the symmetry of the excitation spectrum (and of the eigenstates) upon the change of the sign of $n$: $D_{-n}=D_n$ and $E_{-n}=E_n$, while $\omega_{-n}=\omega_n$. The normalization condition~\eqref{norma2} implies that $|D_n|^2 - |E_n|^2 = 1$.

\section{Density and velocity distribution}\label{sect4}

The number of atoms per unit angle at fixed chemical potential $\mu = g\rho_0$ differs from the expected value $\frac{N_0}{2\pi}$ due to quantum fluctuations and is given by $\nu(\theta) = \langle \hat \psi^\dagger(\theta)\hat \psi(\theta) \rangle$ where the average is taken on the vacuum of the $\hat b_n$ operators, while the total number of particles is $\int_0^{2\pi} d\theta \, \nu(\theta) = N$. The field operator $\hat \psi(\theta)$ has been written as the sum of a classical and a fluctuation part: $\hat\psi(\theta) = \Phi(\theta) +\delta\hat\psi(\theta)$~\eqref{split0}. Expanding to second order in $\delta\hat\psi(\theta)$ and using Eq.~\eqref{expa}, the density profile in the ground state is explicitly obtained:
\begin{equation}
\nu_0(\theta) = |\Phi(\theta)|^2 + \sum_{n\ne 0} \, |E_n(\theta)|^2 = \frac{N_0+ \sum_{n\ne 0} |E_n|^2}{2\pi} =\frac{N}{2\pi}
\end{equation}
which is uniform along the ring.

It is also possible to evaluate the velocity (or angular momentum) distribution in the BEC. The Fourier transform of the field operator reads
\begin{equation}\label{fourier}
\hat \psi_n = \int_0^{2\pi} \frac{d\theta}{\sqrt{2\pi}} \, e^{-i\theta n} \,\hat\psi(\theta) = D_{n}\,\hat b_{n} + E^*_{-n}\,\hat b^\dagger_{-n}
\end{equation}
for $n\ne 0$, while $\hat \psi_0 = \sqrt{N_0}$. Note that, due to the properties of the coefficients $(D_n,E_n)$, $\hat \psi_n$ obeys the canonical commutation relations. This relation can be inverted to give
\begin{equation}\label{bpsi}
\hat b_{n} = D^*_n \hat \psi_n - E^*_n \hat \psi^\dagger_{-n}
\end{equation}
The number of particles with quantized angular momentum $\hbar n \ne 0$ (i.e. with velocity $\frac{\hbar n}{M R_0}$) is given by
\begin{equation}\label{averagen}
G_0(n)=\left\langle \hat \psi^\dagger_n\hat \psi_n \right\rangle = |E_{-n}|^2
\end{equation}

Finally, we can evaluate the correlation function in the ground state
\begin{equation}\label{defgnn}
G_0(n^\prime,n) = \left\langle \hat \psi^\dagger_{n^\prime}\hat \psi_{n^\prime}\hat \psi^\dagger_n\hat \psi_n\right\rangle -
\left\langle \hat \psi^\dagger_{n^\prime}\hat \psi_{n^\prime}\right\rangle\,\left\langle \hat \psi^\dagger_n\hat \psi_n\right\rangle
\end{equation}
By inserting the previous expressions, we find that particles with angular momentum $\hbar n$ and $\hbar n^\prime$ are correlated only if $n^\prime = \pm n \ne 0$:
\begin{equation}\label{gnn}
G_0(n^\prime,n) = \left [\delta_{n,n^\prime}+\delta_{n,-n^\prime}\right] \, E^*_{-n}E_{-n^\prime} D^*_{n}D_{n^\prime}
\end{equation}

It is also possible to obtain the full distribution of angular momenta in the BEC. From Eq.~\eqref{fourier} we can express the ground state of the system, i.e. the vacuum of the $\hat b_n$ ($|0\rangle$), in terms of the vacuum of $\hat \psi_n$ ($|\Psi\rangle$). Starting from Eq.~\eqref{bpsi} we obtain:
\begin{equation}
|0\rangle =
\frac{e^{\,\sum\limits_{n\,>\,0} \frac{E^*_n}{D^*_n} \,\hat \psi^\dagger_n \,\hat\psi^\dagger_{-n}}|\Psi\rangle}{\prod\limits_{n\,>\,0} |D_n|}
\end{equation}
In this form we recognize the squeezed vacuum state, populated by pairs of particles with opposite angular momentum $\pm n$. Expanding the exponential, for each mode $n >0$ the ground state is written as
\begin{equation}
|0\rangle = |D_n|^{-1}\,\sum_{\nu=0}^\infty \left [\frac{E^*_n}{D^*_n}\right]^\nu \,|\nu\rangle
\end{equation}
where $|\nu\rangle$ is the normalized state with $\nu$ particles of angular momentum $n$ and $\nu$ particles of angular momentum $-n$. Therefore the probability $P_n(\nu)$ that the modes ($n,-n)$ are populated by $\nu$ particles each is given by a thermal distribution
\begin{equation}\label{prob}
P_n(\nu) = |D_n|^{-2} \left | \frac{E_n}{D_n}\right |^{2\nu}
\end{equation}

\section{Time dependence: Forced oscillations}\label{sect-time}

Now we study the time evolution of the previously characterized ground state when the ring radius undergoes periodic oscillations:
\begin{equation}\label{rt}
R_t = R_0 \,(1+\epsilon\,\sin\Omega t)
\end{equation}
where $\epsilon \ll 1$ is the amplitude of the perturbation and $\Omega$ its frequency.

The time evolution is investigated adopting the Heisenberg picture, where the state is time independent and coincides with the unperturbed ground state, i.e. the vacuum $|0\rangle$ of the annihilation operators $\hat b_n$~\eqref{inverse}, while the field operators $\hat \psi(\theta,t)$ are now time dependent and satisfy the evolution equation
\begin{equation}
i\hbar\frac{\partial \hat\psi(\theta,t)}{\partial t} = - \left [\hat H,\hat\psi(\theta,t)\right]
\end{equation}
Analogously to the static case, we set
\begin{equation}\label{split}
\hat\psi(\theta,t) = \psi(\theta,t) +\delta\hat\psi(\theta,t)
\end{equation}
where the ``condensate wavefunction'' $\psi(\theta,t)$ is a suitable complex function to be determined. Clearly, the commutation relations are satisfied for any choice of $\psi(\theta,t)$ provided $\delta\hat\psi(\theta,t)$ obeys Eq.~\eqref{commuta} at any time $t$.

Performing the commutations, keeping terms up to first order in $\delta\hat\psi$ and equating the classical and the fluctuation terms separately, we get:
\begin{align}
i\hbar\frac{\partial \psi(\theta,t)}{\partial t} &=
\Big (-\frac{\hbar^2}{2MR_t^2} \frac{d^2}{d\theta^2} + \frac{g}{R_t}\,|\psi(\theta,t)|^2 - g\rho_0\Big)\psi(\theta,t) \label{sch1}\\
i\hbar\frac{\partial \delta\hat\psi(\theta,t)}{\partial t} &= h_2\,\delta\hat\psi(\theta,t) +\frac{g}{R_t}\psi(\theta,t)^2 \delta\hat\psi^\dagger(\theta,t)\qquad\label{sch2}
\end{align}
where terms quadratic in the fluctuations have been neglected, as usual, and
\begin{equation}
h_2= -\frac{\hbar^2}{2MR_t^2} \frac{d^2}{d\theta^2} +2\,\frac{g}{R_t}\,|\psi(\theta,t)|^2 - g\rho_0
\end{equation}
The fluctuation operator can be expressed in terms of the previously defined (time independent) creation and annihilation operators as:
\begin{equation}\label{expand}
\delta\hat\psi(\theta,t) = \sum_n \left \{ D_n(\theta,t)\,\hat b_n + E_n^*(\theta,t)\hat b^\dagger_n\right \}
\end{equation}
where the complex functions $D_n(\theta,t)$ and $E_n(\theta,t)$ obey the usual orthogonality conditions~\eqref{norma} in order to enforce the commutation relations.

The evolution equations of the coefficients are:
\begin{equation}
\begin{split}
i\hbar \frac{\partial D_n(\theta,t)}{\partial t} &= h_2 D_n(\theta,t) + \frac{g}{R_t}\,\psi(\theta,t)^2 E_n(\theta,t) \\
-i\hbar \frac{\partial E_n(\theta,t)}{\partial t} &= h_2 E_n(\theta,t) + \frac{g}{R_t}\,{\psi(\theta,t)^*}^2 D_n(\theta,t)
\end{split}
\end{equation}
to be solved with initial conditions: $\psi(\theta,0) = \Phi(\theta)=\sqrt{\frac{N_0}{2\pi}}$, $D_n(\theta,0) = D_n(\theta)$ and $E_n(\theta,0) = E_n(\theta)$, where $D_n(\theta)$ and $E_n(\theta)$ are defined in Eqs.~\eqref{bogo2}, \eqref{eig}.

Again, the time dependent GPE equation~\eqref{sch1} for the condensate wavefunction is easily solved: The condensate remains uniform at any time with
\begin{equation}
\psi(\theta,t) = \psi(t) = \sqrt{\frac{N_0}{2\pi}} \,e^{i\,\epsilon\,\gamma(t)}
\end{equation}
where $\epsilon\hbar\dot\gamma(t) = g\rho_0 [1-\frac{R_0}{R_t} ]$ which gives $\gamma(t)=\frac{g\,\rho_0}{\hbar\Omega} \,(1-\cos \Omega t)$ for $\epsilon \ll 1$. Also the Bogoliubov modes maintain the same angular dependence:
\begin{align}
\begin{split}
D_n(\theta,t) &= D_n(t)\,\frac{e^{in\theta+i\epsilon\gamma(t)}}{\sqrt{2\pi}} \\
E_n(\theta,t) &= E_n(t)\,\frac{e^{in\theta-i\epsilon\gamma(t)}}{\sqrt{2\pi}}
\end{split}
\end{align}
To first order in $\epsilon$, the equations satisfied by the time dependent amplitudes $D_n(t),E_n(t)$ can be cast in the form:
\begin{equation}\label{dtet}
i\hbar\frac{d}{d t}
\begin{bmatrix}
D_n(t)\\
E_n(t)
\end{bmatrix} =
\begin{bmatrix}
\frac{\hbar^2 n^2}{2MR_0^2} + g\rho_0& g\,\rho_0 \\
-g\,\rho_0& -\frac{\hbar^2 n^2}{2MR_0^2} - g\rho_0
\end{bmatrix}
\begin{bmatrix}
D_n(t)\\
E_n(t)
\end{bmatrix} - \epsilon \,\sin\Omega t \,\clL\,
\begin{bmatrix}
D_n(t)\\
E_n(t)
\end{bmatrix}
\end{equation}
where the perturbation matrix $\clL$ is
\begin{equation}
\clL =
\begin{bmatrix}
\frac{\hbar^2 n^2}{MR_0^2} + g\rho_0& g\,\rho_0 \\
-g\,\rho_0& -\frac{\hbar^2 n^2}{MR_0^2} - g\rho_0
\end{bmatrix}
\end{equation}
Now we expand the unknown functions $D_n(t)$ and $E_n(t)$ on the two right eigenvectors of the matrix~\eqref{eig}. In fact, a basis in the two dimensional space at fixed $n$ requires two linearly independent vectors, which can be chosen as the previously defined positive and the negative frequency right eigenstates:
\begin{equation}\label{expand1}
\begin{bmatrix}
D_n(t)\\
E_n(t)
\end{bmatrix} = \sum_{\eta=\pm} c_n^\eta(t) \,e^{-i\eta\omega_n t}
\begin{bmatrix}
D_n^\eta\\
E_n^\eta
\end{bmatrix}
\end{equation}
where now the stationary eigenvectors, solution of Eq.~\eqref{eig}, also depend on the branch index $\eta$ and are simply related by $D_n^-=E_n^+$, $E_n^-=D_n^+$.
Note that, being the matrix~\eqref{eig} and its eigenvalues real, the eigenvectors $(D_n^\eta,E_n^\eta)$ can always be chosen as real. The normalization condition gives $|c_n^+(t)|^2 - |c_n^-(t)|^2=1$. Substituting into the evolution equation and projecting onto the $\eta^{th}$ eigenmode we get:
\begin{equation}\label{proj}
i\eta\hbar \,\dot c^\eta_n(t) = -\epsilon\sin\Omega t \, \sum_{\eta^\prime=\pm}
\left\{ c^{\eta^\prime}_n(t) \,e^{i\omega_n\,(\eta-\eta^\prime) t}
\, \left [D_n^\eta \,\,-E_n^\eta\right] \,\clL
\begin{bmatrix}
D_n^{\eta^\prime}\\
E_n^{\eta^\prime}
\end{bmatrix}
\right\}
\end{equation}
Note that $\clL$ mixes the positive and negative eigenmode just because the kinetic and the interaction term scale differently with $R_t$. In the quintic GPE both terms are proportional to $(R_0/R_t)^2$ and $\clL$ leaves the two eigenspaces invariant. As previously stressed, this fact is going to kill the Dynamical Casimir Effect in our set-up at strong coupling.

By rescaling time units by $\epsilon$ and defining a ``secular'' time as $\tau = \epsilon\,t$, this equation may be written in the suggestive form
\begin{equation}
\dot c^\eta_n(\tau) = \frac{\eta}{2\hbar}\, \sum_{\eta^\prime=\pm} \left\{ c^{\eta^\prime}_n(\tau)\,
\left [e^{i\frac{\omega_n\,(\eta-\eta^\prime)+\Omega}{\epsilon} \tau} - e^{i\frac{\omega_n\,(\eta-\eta^\prime)-\Omega}{\epsilon} \tau} \right]
\,\left [D_n^\eta \,\,-E_n^\eta\right] \,\clL
\begin{bmatrix}
D_n^{\eta^\prime}\\
E_n^{\eta^\prime}
\end{bmatrix}
\right\}
\end{equation}
showing that rapid oscillations with zero average develop during evolution, unless $\Omega = 2\omega_n$. In this case, the contribution coming from the terms $\eta^\prime=-\eta$ drives the evolution on the ``secular'' time-scale $\tau$. Let us assume that indeed, for a certain mode that we will identify with the label $m$, the resonance condition $\Omega = 2\omega_m$ is satisfied. Keeping only these terms, the secular evolution of the $m^{th}$ mode is then governed by the equations:
\begin{equation}
\dot c^\pm_m(\tau) = K c^{\mp}_m(\tau)
\end{equation}
with a real prefactor $K$ given by
\begin{equation}\label{k}
K = -\frac{\hbar m^2}{2M R_0^2} D_m \,E_m = \frac{\omega_m}{2} \left [1+ \sqrt{1+\left (\frac{\hbar\omega_m}{g\rho_0}\right)^2}\right]^{-1}
\end{equation}
where, again, we dropped the argument $\eta=+$ for the positive branch. The solution, with initial condition $c^+_m(0)=1$ and $c^-_m(0)=0$ is just
\begin{equation}
c^+_m(\tau) = \cosh K\tau \qquad\qquad c^-_m(\tau) = \sinh K\tau
\end{equation}
which satisfies the normalization condition, as requested.

Recalling that the mode $-m$ is degenerate, we find that also this mode grows secularly in time at the same rate. The time dependent coefficients $D_{\pm m}(t)$ and $E_{\pm m}(t)$ are explicitly given by
\begin{equation}\label{dtet2}
\begin{split}
D_{\pm m}(t) &= \cosh K\tau \,e^{-i\omega_m t} D_m + \sinh K\tau \,e^{i\omega_m t} E_m\\
E_{\pm m}(t) &= \cosh K\tau \,e^{-i\omega_m t} E_m + \sinh K\tau \,e^{i\omega_m t} D_m
\end{split}
\end{equation}
while the evolution of the fluctuation operator~\eqref{expand} due to the external forcing grows in time as:
\begin{equation}
\delta\hat\psi(\theta,t) =
\frac{e^{i\epsilon\gamma(t)}}{\sqrt{2\pi}}\,\Big [e^{im\theta} \,D_m(t)\,\hat b_m + e^{-im\theta} \,E^*_m(t)\,\hat b^\dagger_m + e^{-im\theta} \,D_m(t)\,\hat b_{-m} + e^{im\theta} \,E^*_m(t)\,\hat b^\dagger_{-m} \Big]
\end{equation}
The mixing of positive and negative frequency eigenstates is apparent in Eq.~\eqref{dtet2}. Comparing this expansion with the $t=0$ case, we see that the fluctuation operator $\delta\hat\psi(\theta,t)$ is formally given by the same expression of $\delta\hat\psi(\theta,0)$ with the substitution
\begin{equation}\label{modem}
\begin{split}
D_{\pm m} &\to e^{i\gamma(t)}D_{\pm m}(t)\\
E_{\pm m} &\to e^{-i\gamma(t)}E_{\pm m}(t)
\end{split}
\end{equation}
where $D_{\pm m}(t)$ and $E_{\pm m}(t)$ are given by Eq.~\eqref{dtet2}. In the Heisenberg picture the state of the system is unchanged and coincides with the vacuum of all $\hat b_n$ operators, therefore all the results obtained for $t=0$ remain valid at finite $t$ via the substitution~\eqref{modem} for the resonant modes $\pm m$.

The fluctuation operator $\delta\hat\psi(\theta,t)$ can be written in a more compact way by defining two suitable time dependent Bogoliubov modes:
\begin{equation}\label{new}
\hat a^\dagger_{\pm m}(\tau) = \cosh K\tau \,\hat b^\dagger_{\pm m} + \sinh K\tau \,\hat b_{\mp m}
\end{equation}
This is a canonical transformation defining new phonon operators for the two resonant modes. In terms of the new operators the fluctuation operator $\delta\hat\psi(\theta,t)$ has the standard expansion in terms of Bogoliubov eigenmodes, formally identical to the unperturbed case $\epsilon=0$ with the replacements $\hat b_{\pm m} \to \hat a_{\pm m}(\tau)$. The effects of breathing is fully included in the definitions~\eqref{new} of the two new modes. However, the state of the system, being the vacuum of the $\hat b_{\pm m}$ operators is not the vacuum of the new $\hat a_{\pm m}(\tau)$ and therefore in this representation the modes $m$ and $-m$ are populated by phonons. More precisely, the vacuum of the $\hat b_{\pm m}$ operators $|0_b\rangle$ is given by
\begin{equation}
|0_b\rangle = \left [\cosh K\tau\right]^{-1}\,e^{\tanh K\tau \,a^\dagger_m a^\dagger_{-m}} |0_a\rangle
\end{equation}
in terms of the vacuum of the $\hat a_{\pm m}(\tau)$ operators $|0_a\rangle$. The ground state of the model is then a two mode squeezed vacuum state which characterizes several fluctuation related phenomena in different frameworks, like the Hawking effect~\cite{hawking}, inflation and cosmology~\cite{cosmology,inflation1,inflation2} and their ``analogue'' counterparts~\cite{isoard,agullo1,jacobson,agullo2}. The entangled phonons pair have opposite momenta: this is a manifestation of the Dynamical Casimir Effect~\cite{parentani}. Note that, by tracing out one of the two ``twin phonons'', this state gives rise to a thermal density matrix at an effective Casimir temperature $T_C$ which grows on a secular timescale:
\begin{equation}\label{tc}
kT_C = \frac{\hbar\omega_m}{2\log\coth K\tau}
\end{equation}
At large times the effective phonon temperature grows exponentially as $kT_C \sim \frac{\hbar\omega_m}{4} e^{2K\tau}$. In order to estimate $T_C$ we consider a $^{23}$Na condensate of $2.5\,10^5$ atoms in a toroidal trap of average radius $R_0=30\,\mu$m obtained following the procedure described in Ref.~\cite{toro}. The mode $m=1$ is supposed to be excited with oscillations of amplitude $\epsilon=0.1$ and angular frequency $\Omega/2\pi=21.86$ Hz. With these parameters, the growth rate is $K=17.16$ Hz and the effective temperature is about $4$ nK after $1$s but becomes $125$ nK after $2$s.

\section{Density and velocity distribution}\label{sect6}

The rotational symmetry of the ground state is not broken by the time dependent perturbation, therefore the density remains uniform and no undulations are generated by the breathing of the ring.

The average number of particles with given angular momentum $\hbar n$ $G_t(n)=\langle \hat \psi^\dagger_n(t)\hat \psi_n(t) \rangle$ was already evaluated at $t=0$ in Eq.~\eqref{averagen}. At finite $t$ only the modes $\pm m$ are affected by the perturbation and their population is immediately obtained from Eq.~\eqref{averagen} by simple substitution:
\begin{equation}\label{gm}
G_t(\pm m)=|E_m(t)|^2 = \Big | \cosh K\tau \,e^{i\omega_m t} E_m +\sinh K\tau\,e^{-i\omega_m t} D_m \Big |^2
\end{equation}
Both resonant modes show the same population at all times $\langle\nu_+\rangle=\langle\nu_-\rangle=G_t(\pm m)$, which grows on the secular timescale with oscillations at the angular frequency of the breathing $\Omega$. Figure~\ref{figure1} shows the expected fractional occupation of each resonant mode $\frac{G_t(\pm m)}{N_0}$ as a function of time, in the experimental conditions detailed in the previus Section.

\begin{figure}[!htbp]
\includegraphics[width=8cm]{fig1.png}
\caption{Fractional number of particles on each resonant mode as a function of time (in seconds). The chosen parameters correspond to a condensate of $N_0=2.5\,10^5$ sodium atoms in a trap of average radius $R_0=30\,\mu$m. The resonant oscillation frequency is $\Omega/2\pi=2\omega_1/2\pi=21.86$ Hz.}\label{figure1}
\end{figure}

After $1$s, the population of the two modes still represents a small fraction of the total condensate, although a significant increase has occurred. Clearly, the exponential growth cannot last indefinitely because at long times the assumption of weak fluctuations breaks down. Particles that initially are in the state of vanishing angular momentum are set into motion by the resonant drive and populate the levels $\pm m$.

The time evolution of the velocity-velocity correlation function~\eqref{defgnn}, when the system is initially set in its ground state, is easily obtained from Eq.~\eqref{gnn} and reads, for the resonant modes,
\begin{align*}
G_t(\pm m,\pm m) &= |E_m(t)|^2 \,|D_m(t)|^2 \\
&= \left | \cosh K\tau \,e^{-i\omega_m t} D_m + \sinh K\tau \,e^{i\omega_m t} E_m \right |^2
\,\left |\cosh K\tau \,e^{-i\omega_m t} E_m + \sinh K\tau \,e^{i\omega_m t} D_m \right |^2
\end{align*}
A distinctive feature of the squeezed vacuum state is the absence of quantum fluctuations in the unbalance between the populations of the two modes $\Delta\nu = \nu_+ - \nu_-$:
\begin{equation}\label{zero}
\left\langle(\nu_+ - \nu_-)^2\right\rangle = 2\,\left [G_t(m,m)- G_t(m,-m)\right] = 0
\end{equation}
although quantum fluctuations in the population of each mode can be huge:
\begin{equation}\label{variance}
\left\langle\left(\nu_+ - \langle\nu_+\rangle\right)^2\right\rangle = G_t(m,m) = |E_m(t)|^2 \,|D_m(t)|^2
\end{equation}
with $D_m(t)$ and $E_m(t)$ growing in time according to Eq.~\eqref{dtet2}. Therefore, the condition $\langle(\nu_+ - \nu_-)^2\rangle \ll 1$ while $\langle(\nu_+ - \langle\nu_+\rangle)^2\rangle \gg 1$, i.e. $G_t(m,m)- G_t(m,-m) \ll 1$ and $G_t(m,m)\gg 1$, is a hallmark of the analog DCE.


\begin{figure}[!htbp]
\includegraphics[width=7.75cm]{fig2.png}
\caption{Variance of the number of particles on each resonant mode normalized to the square average occupation of the mode as a function of time (in seconds). The chosen parameters correspond to a condensate of $N_0=2.5\,10^5$ sodium atoms in a trap of average radius $R_0=30\,\mu$m. The resonant oscillation frequency is $\Omega/2\pi=2\omega_1/2\pi=21.86$ Hz.}
\label{figure2}
\end{figure}

Figure~\ref{figure2} shows the variance of the number of particles on the resonant mode~\eqref{variance} normalized to the square of the average occupation of the mode $[G_t(m)]^2$. Sharp oscillations of angular frequency $\Omega$ are clearly visible, together with a strong increase on the secular timescale. A quantitative measure of the nature of the squeezed vacuum state is the quantum noise reduction factor, defined as the ratio between the variance of the difference between the populations of the two modes and their average population:
\begin{equation}\label{nrf}
R=\frac{\left\langle (\nu_+-\nu_-)^2\right\rangle - \left\langle (\nu_+-\nu_-)\right\rangle^2}{\langle \nu_+\rangle +\langle\nu_-\rangle}
\end{equation}
In the ideal case represented by the model studied here, this ratio would be zero, but in a real experiment a value less than unity is enough to identify the quantum nature of correlations characterizing the DCE~\cite{graham,nrf}.

The emerging picture would be different if initially the system is set in an excited state of the phonon operators, say due to a small not axisymmetric perturbation. For instance, we can start from a coherent state of phonons of angular momentum $\hbar m$ at $t=0$ breaking rotational symmetry (while modes with $n\ne m$ are in the vacuum state):
\begin{equation}\label{coherent}
|\alpha\rangle = e^{-\frac{|\alpha|^2}{2}} e^{\alpha\hat b_m^\dagger} |0\rangle
\end{equation}
where $\alpha$ is a complex number characterizing the coherent state. The initial average density profile is now non uniform and evolves in time as
\begin{equation}
\left\langle \hat \psi^\dagger(\theta,t)\hat \psi(\theta,t) \right\rangle = \nu_t(\theta) =
\frac{N_0}{2\pi} + \frac{\sqrt{N_0}}{2\pi}\,
\left\{ e^{im\theta} \left[D_m(t)+E_m(t)\right] \,\left\langle \hat b_m\right\rangle + c.c.
\right\}
\end{equation}
to lowest order in the fluctuations. If the resonance condition $\Omega = 2\, \omega_m$ is verified, the coefficients $D_m(t)$ and $E_m(t)$ grow exponentially on the secular timescale leading to
\begin{equation}
\nu_t(\theta) = \frac{N_0}{2\pi} + \frac{\sqrt{N_0}}{\pi}\,\alpha\,\left(D_m+E_m\right)
\,\left[e^{K\tau} \,\cos m\theta\,\cos\omega_m t + e^{-K\tau} \,\sin m\theta\,\sin\omega_m t \right]
\end{equation}
where we set $\alpha$ real. In the absence of external drive, the time evolution would be given by the same expression with $K=0$, i.e. a standing wave $\cos(m\theta-\omega_m t)$ on top of a uniform condensate. This result shows that while the solution of the unperturbed GPE is dynamically stable according to the standard criteria~\cite{stab}, quantum fluctuations give rise to an instability on the secular timescale: if initially the state has a small undulation of the correct periodicity, such a modulation is exponentially amplified by quantum fluctuations. This is the familiar mechanism of parametric resonance~\cite{nori}, present also in a classical framework. As an example, in Figure~\ref{figure0} we show the density profile after $1,2$ and $3$ seconds
in a system defined by the same parameters previously quoted. The initial density modulation is set at $1\%$. After $3$s the amplitude has become
as large as $50\%$.

\begin{figure}[!htbp]
\includegraphics[width=8cm]{fig0.png}
\caption{Density profile $\nu_t(\theta)$, normalized to the initial value $\nu_0(\theta)$ for a condensate of $N_0=2.5\,10^5$ sodium atoms in a trap of average radius $R_0=30\,\mu$m. At $t=0$ the system is set in a coherent
state characterized by a small density modulation of $1\%$.
(black curve). Modulations of larger amplitudes refer to increasing evolution times: $1$s (green), $2$s (blue) and $3$s (red).}
\label{figure0}
\end{figure}

Remarkably, also in this case, the noise reduction factor is always smaller than unity and decreases due to the amplification of the fluctuations signaling the quantum nature of the stimulated dynamical Casimir effect:
\begin{equation}
R_s = \left [1 + 2\,|E_m(t)|^2 \,\left (1 + |\alpha|^{-2}\right)\right]^{-1}
\end{equation}
where $E_m(t)$ is given in Eq.~\eqref{dtet2}. This result formally coincides with the one derived in the context of quantum optics~\cite{bondani3,bondani}.

In a slightly different set-up, both resonant modes $\pm m$ are initially excited and a stationary wave is generated in the condensate before switching-on the periodic modulation. Formally this corresponds to set both modes in a coherent state~\eqref{coherent} and let the system evolve. Again, if the resonance condition is satisfied, the amplitude of the standing wave increases exponentially on the secular time-scale $\tau=\epsilon \,t$. Defining $\delta \nu_t(\theta) = \nu_t(\theta) - \frac{N_0}{2\pi}$ we find
\begin{equation}\label{deltanu}
\delta \nu_t(\theta) = 2\frac{\sqrt{N_0}}{\pi}\,\alpha\,(D_m+E_m)
\,e^{K\tau} \cos m\theta\,\cos\omega_m t
\end{equation}
for real $\alpha$. Also the population of each resonant mode $\langle \nu_{\pm}\rangle=G(\pm m)$ grows exponentially, while the variance of the population difference $\langle (\nu_+ -\nu_-)^2\rangle$ remains finite and time independent, equal to $2\,|\alpha|^2$. However, in this case, the noise reduction factor~\eqref{nrf} starts with a value larger than unity while crosses-over to $R<1$ due to parametric amplification. This result means that the criterion for detecting quantum correlations is reached only after a sufficiently long pumping time.

Not surprisingly, finite temperature effects are able to hide the quantum correlations, at least in the first stages of the amplification. In the absence of external perturbations, the noise reduction factor is given by:
\begin{equation}\label{rt}
R_T = \left [1 - e^{-\beta\hbar\omega_m} + 2\,|E_m(t)|^2 \,\sinh (\beta\hbar\omega_m)\right]^{-1}
\end{equation}
where $\beta = 1/k_BT$ and again $E_m(t)$ is given in Eq.~\eqref{dtet2}. This expression agrees with the analogous result in the quantum optics framework~\cite{bondani3,bondani2}.
A plot of the predicted time dependence of $R_T$ for three temperatures in the experimental set up previously considered is shown in Figure~\ref{figure3}. Strong oscillations are present at short times but after $2$s the noise reduction factor is smaller than unity for most of the time. Interestingly, temperature mainly affects the height of the sharp peaks but does not severely limits the detection of the spontaneous DCE in this set-up.

\begin{figure}[!htbp]
\includegraphics[width=8cm]{fig3.png}
\vspace*{-15pt}
\caption{Noise reduction factor $R_T$ (see Eq.~\eqref{rt}) as a function of time (in seconds). The chosen parameters correspond to a condensate of $N_0=2.5\,10^5$ sodium atoms in a trap of average radius $R_0=30\,\mu$m. The resonant oscillation frequency is $\Omega/2\pi=2\omega_1/2\pi=21.86$ Hz. The main panel shows the results for three temperatures: $T=50$ nK (black), $T=100$ nK (blue) and $T=200$ nK (red). In the inset, the data at $T=50$ nK are shown in an extended temporal range.}
\label{figure3}
\end{figure}


\section{The classical counterpart}\label{sect7}

Here we describe the time evolution of the same set-up in terms of a purely classical model, showing that very similar results are obtained as far as the parametric amplification is concerned. However, crucial differences emerge in the velocity distribution function, which represents the hallmark of quantum fluctuations.

The GPE is known to be equivalent to the Euler equation for a classical perfect fluid~\cite{stringa} with the inclusion of an additional ``quantum pressure term''. In our case, Eq.~\eqref{gpe1d} can therefore be written as the continuity and Euler equation for the number density $\rho(\theta,t)$ and fluid velocity $v(\theta,t)$:
\begin{align}
\frac{\partial (\rho R_t)}{\partial t } + \frac{\partial (\rho v)}{\partial \theta } &= 0\label{eqc} \\
\rho R_t \frac{\partial v}{\partial t } + \rho v \frac{d R_t}{d t } + \rho v \frac{\partial v}{\partial \theta } &= -M^{-1}\,\frac{\partial p }{\partial \theta }\label{euler}
\end{align}
where the quantum pressure contribution has been disregarded and the equation of state of the fluid is of the polytropic form $p(\rho) = C \,\rho^\gamma$. The GPE equation is recovered with the choice $C=\frac{g}{2}$ and $\gamma=2$. The second contribution is the ``Hubble friction term''~\cite{jacobson} which can be eliminated by introducing the number of particles per unit angle $\nu_t(\theta) = \rho(\theta,t) \,R_t$ and the rescaled velocity $u(\theta,t)=v(\theta,t)\,R_t$ in place of the particle density $\rho$ and velocity $v$. The pair of equations~\eqref{eqc}, \eqref{euler} describe a classical, ideal, polytropic fluid in a ring whose radius $R_t$ periodically oscillates according to Eq~\eqref{rt}. In the static case $\epsilon=0$, the trivial solution $\nu^0=\frac{N_0}{2\pi}$, $u=0$ represents the phonon vacuum state. Upon turning on the oscillations, the density is perturbed by $\nu_t(\theta) = \nu^0 + \delta \nu_t(\theta)$, while a non-vanishing fluid velocity $u(\theta,t)$ sets in. Linearizing the equations in $\delta \nu$ and $u$ and setting as initial condition a weak standing wave of the form~\eqref{deltanu} $\delta \nu_0(\theta) = \delta \nu_0 \,\cos(\theta m)$ the solution reads
\begin{equation}\label{perturba}
\delta \nu_t(\theta) = \delta \nu_0 \,e^{K\tau}\cos(\theta m) \,\cos(\omega_m t)
\end{equation}
where the eigenfrequency $\omega_m$ is given by
\begin{equation}
\omega_m^2 = \frac{C\,\gamma\,\left(\nu^0\right)^{\gamma-1}}{R_0^{\gamma+1} M} \, m^2
\end{equation}
and coincides with the long wavelength limit of the Bogoliubov expression~\eqref{omega}. In obtaining~\eqref{perturba}, the usual resonance condition $\Omega = 2\omega_m$ has been assumed. Interestingly, the parametric amplification is governed by the rate
\begin{equation}
K = \frac{3-\gamma}{4} \,\omega_m
\end{equation}
which is positive for $\gamma < 3$ and agrees with the result of the GPE analysis~\eqref{k} for $\gamma=2$, but vanishes for $\gamma=3$, which indeed corresponds to the Tonks--Girardeau strong coupling limit. Therefore, the classical result for the evolution of the density profile~\eqref{perturba} coincides with the quantum calculation~\eqref{deltanu} showing that the phenomenon of parametric amplification does not reveal the quantum features of the system. Instead, the growing peaks in the particle velocity distribution $G_t(n)$~\eqref{gm} at the resonant modes $n=\pm m$ is a purely quantum feature which has no analog in a classical continuum theory.

\section{Summary and Conclusions}\label{sect8}

The analysis previously carried out proves that the ``breathing ring geometry'' is a promising set-up for an unambiguous detection of the analog DCE in BEC. While \emph{stimulated} DCE was indeed observed in the experiment of Ref.~\cite{jaskula}, the ring geometry will allow to detect also the \emph{spontaneous} effect, implying the amplification of vacuum fluctuations. Contrary to previous proposals~\cite{jaskula,demler} the imposed periodic perturbation does not directly couple to the condensate giving rise to mechanical waves propagating through the fluid. When the resonance condition $\Omega = 2\omega_m$ is satisfied, the only effect of the oscillations in the radius of the torus is the creation of entangled pairs of phonons of opposite momenta. Thereby, the population of the resonant modes grows exponentially in time, while the condensate, initially in the unperturbed ground state, is not set into motion neither displays oscillations: Only the velocity distribution of the atoms in the condensate is affected, thereby providing a clear-cut demonstration of the amplification of vacuum fluctuations predicted by the DCE. The quantum nature of this phenomenon is proved by measuring the variance of the number of phonons on the two modes through the noise reduction factor.

The analysis of this phenomenon is carried out truncating the fluctuation hamiltonian to second order in the Bogoliubov modes. Therefore, coupling between modes is neglected. As a consequence, no damping mechanism is present and quantum excitations have infinite lifetime. To include finite lifetime effects, notably the expected Landau damping of phonons, cubic terms should be included perturbatively in the Bogoliubov hamiltonian~\cite{stringari}. Back-reaction effects are expected to become relevant in the late stages of parametric amplification~\cite{butera}, when the perturbation induced by DCE starts to affect the background vacuum state. The inclusion of weak dissipation mechanisms in the DCE was studied in Ref.~\cite{parentani}.

The careful determination of the velocity correlation function can prove the occurrence of the analog DCE via the suppression of fluctuations in the difference between the populations of the two resonant modes. Conversely, if a density modulation of the right wavelength is initially present, the breathing amplifies such an undulation on a secular timescale leading to clearly visible effects on the density profile of the condensate due to the classical phenomenon of parametric resonance. This experiment may be modeled by starting with a coherent, rather than vacuum, state: also in this case parametric amplification progressively lowers the noise reduction factor, leading to a quantitative evidence of entanglement in the phonon state. Although finite temperature effects may hide the quantum nature of correlations, parametric amplification eventually prevails and the noise reduction factor crosses unity at late times.

The phonon dynamics in ring geometry has been also investigated in a different context in Ref.~\cite{jacobson}, where the case of a rapidly expanding condensate is studied and the possible emergence of topological excitations is considered. The evolution equation for the phonon amplitudes derived in that paper coincides with an approximation of the time dependent Bogoliubov equations obtained in Section~\ref{sect-time}.

As an alternative set-up, the ring-shaped condensate may be periodically triggered by modulating the frequency of the transverse trap~\cite{butera}. Also in this case, the symmetry is preserved and we expect that the same analysis previously performed for the breathing ring applies.

\section*{Acknowledgements}
The present investigation was stimulated and motivated by Renaud Parentani during his visits to Como. In particular, the adoption of a toroidal geometry for the study of the Dynamical Casimir Effect was a result of the enlightening discussions we had with Renaud and of his insightful comments. It is a pleasure also to thank Maria Bondani for thoughtful discussions about the non-classicity criteria. This work benefited from scientific exchanges during the workshop ``Quantum simulation of gravitational problems on condensed matter analog models'', organized by ECT* (European Centre for Theoretical Studies in Nuclear Physics and related areas) in Trento (June 19-23, 2023).

%The authors do not work for, advise, own shares in, or receive funds from any organisation that could benefit from this article, and have declared no affiliation other than their research organisations.


%\nocite{*}

\section*{Declaration of interests}
The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

\printbibliography
\end{document}
