${\LARGE MANY-BODY\,EFFECTS\,AND\,QUANTUM\,\ FLUCTUATIONS\vskip 12.0pt plus
4.0pt minus 4.0pt}$

${\LARGE\ FOR\ DISCRETE\,\,TIME\,CRYSTALS\vskip 12.0pt plus 4.0pt minus
4.0pt}$

$\,{\LARGE IN\,}\ {\LARGE BOSE-EINSTEIN\,CONDENSATES}$

Jia Wang1, Peter Hannaford2 and Bryan J. Dalton1

1Centre for Quantum Science and Technology Theory, Swinburne University of
Technology, Melbourne, Victoria 3122, Australia

2Optical Sciences Centre, Swinburne University of Technology, Melbourne,
Victoria 3122, Australia

Abstract We present a fully comprehensive multi-mode quantum treatment based
on the truncated Wigner approximation (TWA) to study many-body effects and
effects of quantum fluctuations on the formation of a discrete time crystal
(DTC) in a Bose-Einstein condensate (BEC) bouncing resonantly on a
periodically driven atom mirror. Zero-range contact interactions between the
bosonic atoms are assumed. Our theoretical approach avoids the restrictions
both of mean-field theory, where all bosons are assumed to remain in a single
mode, and of time-dependent Bogoliubov theory, which assumes boson depletion
from the condensate mode is small. We show that the mean-field and time-
dependent Bogoliubov approaches can be derived as approximations to the TWA
treatment. Differing initial conditions, such as a finite temperature BEC, can
also be treated. For realistic initial conditions corresponding to a harmonic
trap condensate mode function, our TWA calculations performed for period-
doubling agree broadly with recent mean-field calculations for times out to at
least $2000$ mirror oscillations, except at interaction strengths very close
to the threshold value for DTC formation where the position probability
density differs significantly from that determined from mean-field theory. For
typical attractive interaction strengths above the threshold value for DTC
formation and for the chosen trap and driving parameters, the TWA calculations
indicate a quantum depletion due to quantum many-body fluctuations of less
than about two atoms out of a total of $600$ atoms at times corresponding to
$2000$ mirror oscillations, in agreement with time-dependent Bogoliubov theory
calculations. On the other hand, for interaction strengths very close to the
threshold value for DTC formation, the TWA calculations predict a large
quantum depletion - as high as about $260$ atoms out of $600$. We also show
that the mean energy per particle of the DTC does not increase significantly
for times out to at least $2000$ mirror oscillations and typically oscillates
around an average value close to its initial value; so TWA theory predicts the
absence of thermalisation. Finally, we find that the dynamical behaviour of
our system is largely independent of whether the boson-boson interaction is
attractive or repulsive, and that it is possible to create a stable DTC based
on repulsive interactions.

##  1 Introduction

Discrete time crystals (DTC) are periodically driven non-equilibrium states of
quantum many-body systems that spontaneously break discrete time-translation
symmetry due to particle interactions and start to evolve with a period $s$
times $(s=2,3,..)$ longer than the period $T\ $of the drive [1], [2], [3],
[4], [5]. Such time crystals are predicted to be robust against external
perturbations and to persist indefinitely in the thermodynamic limit of large
$N$, in analogy with space crystals. A number of different platforms for
creating discrete time crystals have been proposed [1], [2], [3], [4].
Preliminary experimental evidence for the realization of DTCs has been
reported for a number of platforms, including a spin chain of interacting ions
[6], nitrogen-vacancy spin impurities in diamond [7], nuclear spins in organic
molecules [8] and ordered ADP crystals [9], [10], and a superfluid Bose-
Einstein condensate of ultra-cold atoms [11], [12], [13]. Overviews of the
topic of time crystals are given in [14], [15], [16], [5].

Another platform that has been proposed for a discrete time crystal and which
is yet to be realized experimentally involves the use of a Bose Einstein
condensate (BEC) of weakly interacting bosonic atoms bouncing resonantly on a
periodically driven mirror [1], [17], [18], [19]. A many-body mean-field
approach with the condensate wave-function expanded as a linear combination of
$s$ Wannier functions (which are defined in the time domain) has been employed
to study such discrete time crystals in the case of period-doubling $(s=2)$
[1], [19], period-quadrupling $(s=4)$ [19] and higher periodicities $s=10-100$
[17], [18]. Also, a time-dependent Bogoliubov theory has been used to treat
the $s=2$ case [19] .

A concern in applying a mean-field (single-mode) or a few-mode approach to
study time crystals and discrete time-translation symmetry breaking is whether
the lack of thermalisation and decay of the condensate in such studies is an
artefact imposed by the adopted approximations [16]. Thermalisation refers to
a common phenomenon that many-body systems subjected to periodic driving
eventually reach infinite temperature, where no time crystal can exist. On the
other hand, describing a many-body system with no more than a few modes
precludes thermalisation and quantum depletion of the condensate by
assumption, raising concerns that the predicted time crystal may not exist
once one considers a more rigorous multi-mode approach.

Features that a theory of DTC should take into account include: the system
being many body, the possible effects of there being many modes that the
bosons could occupy and whether one or two modes could dominate, the roles of
both interactions and driving, the evolution dynamics being quantum and
allowing for a steady periodic non-equilibrium state to appear with period a
multiple of the drive period, whether or not thermalisation occurs, and the
sensitivity of DTC behaviour to changes in the initial conditions - including
the effect of the initial BEC temperature. The truncated Wigner approximation
(TWA) can take into account all of these features - whereas (for example) a
mean-field theory only allows for the bosons to occupy one mode, and time-
dependent Bogoliubov theory does not allow for substantial depletion of the
condensate mode.

In this paper we present a full theoretical many-body study of conditions for
creating discrete time crystals in a Bose-Einstein condensate which is allowed
to bounce resonantly on an oscillating mirror in a gravitational field as
proposed in [1]. We use the TWA to predict/ confirm the presence of the non-
trivial quantum phenomena of DTC using a more advanced and accurate method
than in previous studies [1], [19]. The main preparation process considered is
that of a weakly interacting BEC in a harmonic trap, strongly confined in the
transverse direction. We study many-body effects and effects of quantum
fluctuations on the evolution for the case of period-doubling $(s=2)$ and for
a zero temperature BEC. Processes that could prevent a DTC being created, such
as depletion of bosonic atoms from the condensate mode, are allowed for, and
the TWA approach is used to examine whether thermalisation - which also
destroys a DTC - occurs in experimentally accessable time periods. We compare
our results with previous studies [1], [19] where the theoretical approach is
based on mean-field theory and the Gross-Pitaevskii equation (GPE), and hence
assumes that the bosonic atoms remain in a single condensate mode, and on
time-dependent Bogoliubov theory, where the depletion from the condensate mode
is assumed to be small.

The theoretical framework outlined in the present paper is intended to provide
the basis for more extensive future studies, including for other periodicities
$sT$ and for non-zero temperature. We also wish to explore more fully the
conditions for $sT$ periodicity than is possible in this initial paper, such
as choices of the boson-boson interaction strength, different initial
conditions, the stability of the time crystal, etc, but using a theory which
is not restricted to assuming all bosons remain in a single condensate mode
(such as in mean-field theory) or that the quantum fluctuations from the
condensate field are small (as in time-dependent Bogoliubov theory). For these
reasons, and because other expositions of some of the theoretical background
do not cover in one paper all the background material, we have set out a full
description of our approach - including the description of the prepared
initial state via Bogoliubov theory as well as the evolution of the bouncing,
driven BEC. Our approach assumes that the total boson number $N$ is large, as
is usual for a BEC.

The truncated Wigner approximation (TWA) [20], [21], [22], [23], [24], is a
phase-space theory method which is a now well-established approach for
treating dynamical behaviour in cold bosonic gases. We expect it to be more
reliable than either the mean-field theory or time-dependent Bogoliubov
theory, since these approaches can be derived from TWA as approximations (see
Appendix 6). The TWA involves both a field-theory approach, where the BEC is
treated as a single quantum field, and a mode-theory approach, where the
separate quantum modes or single-particle states that the modes could occupy,
are used. The modes that describe the evolution of the BEC are either Floquet
modes, which describe single particles moving in a time-periodic potential, or
gravitational modes, which describe single particles moving in a static
gravitational potential. Both in the single quantum field and the separate
mode approaches, observable quantities at each time can be expressed as phase-
space averages involving a Wigner distribution functional (function) or as a
stochastic average involving stochastic fields or mode amplitudes. For phase-
space averages the observable quantities are first expressed in terms of
symmetrically ordered combinations of the quantum operators (see Ref. [22]).
Similarly, evolution equations for the quantum state can be described by
Fokker-Planck equations for the Wigner functional (function), or as in the
present paper by equivalent Ito stochastic field equations in the field-theory
version or Ito stochastic differential equations for stochastic phase-space
variables in the mode-theory version [21], [22]. This equivalence ensures that
phase space and stochastic averages for observable quantities are equal. The
Ito equations are deterministic in the present case, but depend on initial
values for the phase-space variables - which are specified by a stochastic
distribution determined from the initial state of the BEC just after it is
prepared. The validity of the TWA is discussed in Refs. [21], [23], as well as
in Section 2.2.4 and Appendix 7.

The preparation of the initial state of the BEC is described by time-
independent Bogoliubov theory [22], [24], [25], [26], [27] in which the BEC is
formed in a harmonic trap potential and the condensate wave function is
obtained from a time-independent Gross-Pitaevskii equation (GPE). The
condensate is assumed to be in a Fock state. The non-condensate field is
required to be orthogonal to the condensate mode, and non-condensate modes are
described both via an arbitrary set of mode functions chosen to be orthogonal
to the condensate mode and also via Bogoliubov mode functions determined from
Bogoliubov-de Gennes equations [25]. There is a pair of Bogoliubov mode
functions for each mode frequency, which satisfy biorthogonality conditions.
Separate Wigner distribution functions for the Bogoliubov modes and for the
condensate mode are obtained when the density operator for the initial state
factorises into condensate and non-condensate density operators \- as in the
Bogoliubov approximation. Stochastic phase-space variables are introduced for
both the condensate mode and the Bogoliubov modes, and the stochastic
properties of these phase-space variables are determined for the case where
the Bogoliubov modes are unoccupied and all bosons are in the condensate mode
in a phase-invariant state - such as would occur at zero temperature. The
initial stochastic phase-space variables for the Floquet or gravitational
modes (and their statistical properties) are determined from stochastic phase-
space variables for the condensate and Bogoliubov modes via matching the
stochastic field functions expanded either in terms of Floquet (or
gravitational) modes or in terms of the condensate and (non-condensate)
Bogoliubov modes. To help understand the dynamical processes involved, we also
apply the TWA treatment to a different initial condition, where the condensate
mode function is described as a linear combination of Wannier-like states.
Here, the non-condensate modes are all orthogonal to the condensate mode, but
not chosen as Bogoliubov modes.

Expressions for various quantities of physical interest are obtained as a
function of evolution time, including the position probability density (PPD);
the spatial position quantum correlation function (QCF) - which describes the
spatial coherence of the BEC; the mean values of one-body projectors (OBP) -
which describe the probabilities of finding a boson in the original condensate
mode or Wannier modes; and the quantum depletion - which describes the loss of
bosons from the condensate mode. The TWA treatment is compared with that based
entirely on mean-field theory as treated via the time-dependent version of the
Gross-Pitaevskii equation. Numerical calculations are performed for evolution
times out to $2000$ mirror oscillations and for a wide range of interaction
strengths as well as the two different initial conditions.

We find that for attractive interactions the combination of boson-boson
interaction and periodic driving enables a stable discrete time crystal to be
created for a realistic harmonic trap initial condition when the magnitude of
the interaction parameter exceeds a _threshold_ value, $|gN|>0.012$ for our
choice of parameters (see Table 1). In the absence of driving, we find that
the interaction strength corresponding to that required to create a DTC is
strong enough to couple many modes, which is indicated by a substantial
quantum depletion. However, with driving, the quantum depletion is strongly
suppressed and the thermalisation is quenched, implying that the absence of
our system’s thermalisation is a genuine many-body effect and a direct
consequence of driving in the presence of a sufficiently strong interaction.
We also find that the dynamical behaviour of our system is largely independent
of whether the boson-boson interaction is attractive or repulsive.

We also find that our TWA results agree well with the mean-field GPE
calculations except at interaction strengths very close to the threshold value
for creation of a single stable wave-packet and formation of a discrete time
crystal. We show in Appendix 6 that the mean-field and time-dependent
Bogoliubov approaches are derivable as approximations to the TWA treatment.
For realistic initial conditions corresponding to the harmonic trap condensate
mode function and typical interaction strengths required to produce a DTC, our
TWA calculations indicate a quantum depletion less than about two atoms out of
a total of $600$ atoms at times corresponding to $2000$ mirror oscillations,
in agreement with recent many-body calculations based on a time-dependent
Bogoliubov approach [19]. However, for interaction strengths very close to the
threshold value for DTC formation, the quantum depletion due to the quantum
fluctuations can be as high as about $260$ atoms out of a total of $600$ atoms
at times corresponding to $2000$ mirror oscillations. Here, the time-dependent
Bogoliubov approach would break down, since Bogoliubov theory assumes that the
non-condensate field is relatively small. In the absence of driving and for a
relatively strong interaction ($0.05\leq|gN|\leq 0.1$), we find that the
quantum depletion can be more than $220$ atoms out of a total of $600$ atoms
at times corresponding to $2000$ mirror oscillations, indicating that driving
in the presence of a sufficiently strong interaction needed to form a DTC
supresses the quantum depletion associated with the quantum many-body
fluctuations.

In Section 2, we describe the phase-space theory used to study the evolution
of the periodically driven many-body system, which includes details of the
truncated Wigner approximation (Sect. 2.2), the periodically driven BEC with
stochastic phase amplitudes for separate modes (Sect. 2.3), the preparation of
a BEC with condensate and Bogoliubov modes (Sect. 2.4), and the quantum state
for the initial BEC (Sect. 2.6). In Section 3, we present numerical results of
the TWA calculations; in Section 4 we discuss the significance of our results
in comparison with previous calculations; and in Section 5 we present our
conclusions and suggest future directions. Details of the theory are set out
in Appendices, which are included as Online Supplementary Material.

##  2 Phase Space Theory for Evolution of a Periodically Driven BEC

The effect of quantum fluctuations due to many-body effects on the predictions
of time crystal behaviour in a BEC system may be treated using phase-space
methods. If the truncated Wigner $W^{+}$ distribution functional [20], [21],
[24], [22] is used, the resulting Ito stochastic field equations are similar
to the Gross-Pitaevski equation that arises using the mean-field
approximation, but now allow for quantum fluctuations due to the fields being
stochastic rather than deterministic. The treatment outlined below is based on
the $1D$ approximation, as in the papers of Sacha et al. [1], [17], [18],
[19].

###  2.1 Hamiltonian

The _Hamiltonian_ for the $1D$ model for the periodically driven BEC is given
in terms of the field operators $\hat{\Psi}(z)$, $\hat{\Psi}(z)^{{\dagger}}$
for the annihilation, creation of a bosonic atom of mass $m$ at position $z$
as [17]

$\displaystyle\hat{H}\mbox{\rule{-1.42262pt}{0.0pt}}$
$\displaystyle=\mbox{\rule{-2.84526pt}{0.0pt}}\int\mbox{\rule{-2.84526pt}{0.0pt}}dz\left(\frac{\hbar^{2}}{2m}\frac{\partial}{\partial
z}\hat{\Psi}(z)^{{\dagger}}\mbox{\rule{-1.42262pt}{0.0pt}}\,\mbox{\rule{-1.42262pt}{0.0pt}}\frac{\partial}{\partial
z}\hat{\Psi}(z)+\hat{\Psi}(z)^{{\dagger}}V(z,t)\hat{\Psi}(z)\right.$
$\displaystyle\hskip
56.9055pt\left.+\frac{g}{2}\hat{\Psi}(z)^{{\dagger}}\hat{\Psi}(z)^{{\dagger}}\hat{\Psi}(z)\hat{\Psi}(z)\right)$
(1) $\displaystyle=\widehat{{\small K}}+\widehat{{\small V}}+\widehat{{\small
U}}$ (2)

and is the sum of kinetic energy, time-dependent potential energy and
interaction energy terms.

The potential energy contains a driving term at frequency $\omega=2\pi/T$ and
is given by

$V(z,t)=mg_{E}\,z(1-\lambda\cos\omega t)$ (3)

This corresponds to an initially prepared BEC being allowed to fall in a
vertical gravitational field ($g_{E}$ is the gravitational acceleration) and
allowed to bounce off an atom mirror oscillating with a period $T$ and
amplitude $\lambda$. The system is described in the oscillating frame.

The interaction energy is based on a zero-range contact interaction between
the bosonic atoms. The $1D$ coupling constant is
$g=2\hslash\omega_{\bot}a_{s}$ \- which is proportional to the s-wave
scattering length $a_{s}$ and the oscillation frequency $\omega_{\bot}$ for
the BEC atoms in a transverse trap - may be tuned via a Feshbach resonance.

###  2.2 Truncated Wigner Approximation

In this section we outline the _truncated Wigner approximation_ (TWA) for
phase-space field theory and set out basic expressions for some quantities of
physical interest such as the position probability density and the first order
quantum correlation function.

####  2.2.1 Functional Fokker-Planck Equation - W+ Distribution

The _functional Fokker-Planck equation_ (FFPE) for the $W^{+}$ _Wigner
distribution functional_ $W[\mathbf{\psi}(z)]$ which represents the quantum
density operator $\widehat{\rho}$ can be obtained via applying the
correspondence rules in conjunction with the Liouville-von Neumann equation
for the density operator. The contributions can be written as a sum of terms
from $\widehat{{\small K}},\widehat{{\small V}},\widehat{{\small U}}$. A
derivation for $3D$ is presented in Ref [22] (see Section 15.1.6 and Appendix
I). The field operators $\hat{\Psi}(z)$, $\hat{\Psi}(z)^{{\dagger}}$ are
represented by two unrelated time-independent c-number fields $\psi(z)$ and
$\psi^{+}(z)$, respectively. For short,
$\mathbf{\psi}(z)\equiv\psi(z),\psi^{+}(z)$, and $W[\mathbf{\psi}(z)]$ is time
dependent.

The kinetic energy term is

$\displaystyle\left(\frac{\partial}{\partial t}W[\mathbf{\psi}(z)]\right)_{K}$
$\displaystyle=\frac{i}{\hslash}\int
dz\,\left[\left\\{\frac{\delta}{\delta\psi(z)}\left(-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial
z^{2}}{\small\psi(z)\,}W[\ \mathbf{\psi}(r)]\right)\right\\}\right.$
$\displaystyle\left.+\left\\{\frac{\delta}{\delta\psi^{+}(z)}\left(\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial
z^{2}}{\small\psi}_{u}^{+}{\small(z)\,}W[\
\mathbf{\psi}(r)]\right)\right\\}\right]$ (4)

and only contributes to the drift term in the FFPE.

The potential energy term is

$\displaystyle\left(\frac{\partial}{\partial t}W[\mathbf{\psi}(z)]\right)_{V}$
$\displaystyle=\frac{i}{\hslash}\int{\small
d}z\,\left[\left\\{\frac{\delta}{\delta\psi(z)}\left((V(z,t))\psi(z)\,W[\
\mathbf{\psi}(r)]\right)\right\\}-\left\\{\frac{\delta}{\delta\psi^{+}(z)}\left((V(z,t))\psi^{+}(z)\,W[\
\mathbf{\psi}(z)]\right)\right\\}\right]$

and also only contributes to the drift term in the FFPE.

The boson-boson interaction term is

$\displaystyle\left(\frac{\partial}{\partial t}W[\mathbf{\psi}(z)]\right)_{U}$
$\displaystyle{\small=}\frac{i}{\hslash}g\int{\small
d}z\,\left[\left\\{\frac{\delta}{\delta\psi(z)}(\psi^{+}(z)\psi(z)-\delta_{C}(z,z))\,\psi(z)\,W[\mathbf{\psi}(z)]\right\\}\right.$
$\displaystyle\left.-\left\\{\frac{\delta}{\delta\psi^{+}(z)}(\psi^{+}(z)\psi(z)-\delta_{C}(z,z))\,\psi^{+}(z)\,W[\mathbf{\psi}(z)]\right\\}\right.$
$\displaystyle\left.-\frac{1}{4}\left\\{\frac{\delta}{\delta\psi(z)}\frac{\delta}{\delta\psi(z)}\frac{\delta}{\delta\psi^{+}(z)}\,\psi(z)\,W[\mathbf{\psi}(z)]\right\\}\right.$
$\displaystyle\left.+\frac{1}{4}\left\\{\frac{\delta}{\delta\psi^{+}(z)}\frac{\delta}{\delta\psi^{+}(z)}\frac{\delta}{\delta\psi(z)}\,\psi^{+}(z)\,W[\mathbf{\psi}(z)]\right\\}\right]$
(6)

where for two position coordinates $z$ and $z^{\\#}$

$\delta_{C}(z,z^{\\#})={\displaystyle\sum\limits_{k}}\phi_{k}(z)\phi_{k}(z^{\\#})^{\ast}$
(7)

with the $\phi_{k}(z)$ being any set of suitable orthonormal mode functions
satisfying the condition
${\displaystyle\int}dz\,\phi_{k}(z)^{\ast}\phi_{l}(z)=\delta_{k,l}$. A unitary
change in these modes does not change $\delta_{C}(z,z^{\\#})$.

The exact FFPE involves only first and third-order functional derivatives -
there are no second-order diffusion terms. In the truncated Wigner $W^{+}$
treatment, the third-order derivatives are discarded so only the first-order
drift terms remain. This proceedure is justisfied on the basis that the region
of phase space that is most important is where
$\psi(z),\psi^{+}(z)\sim\sqrt{N}$, where $N$ is the number of bosons present.
As applying a further derivative gives a contribution of order $1/\sqrt{N}$
smaller than the previous term, discarding the third-order derivative terms is
justified when $N$$\gg 1$. The validity of the TWA is discussed in Refs. [21],
[23], as well as in Section 2.2.4 and Appendix 7.

####  2.2.2 Ito Stochastic Field Equations for Time Evolution

The Ito stochastic field equations (SFE) Ito equations are equivalent to the
functional Fokker-Planck equation for the distribution functional assuming the
latter only involve up to second-order functional derivatives. There are no
further approximations involved. This equivalence is demonstrated in Chapter
14 of Ref [22] and is based on the requirement that the phase space average of
any function of the phase space field functions obtained from the Wigner
distribution functional is the same as the stochastic average of the same
function of stochastic field functions. The form of the Ito stochastic
equations that enables this equivalence is obtained from the terms in the
FFPE.

The Ito (SFE) corresponding to the functional Fokker-Planck equation for the
$W^{+}$ distribution functional are as follows: (see Section 15.1.8 in Ref.
[22])

$\displaystyle\frac{\partial}{\partial t}\widetilde{\psi}(z,t)$
$\displaystyle=-\frac{i}{\hslash}\left[-\frac{\hslash^{2}}{2m}\frac{\partial^{2}}{\partial
z^{2}}\widetilde{\psi}(z,t)+V(z,t)\widetilde{\psi}(z,t)+g\\{\widetilde{\psi}^{+}(z,t)\widetilde{\psi}(z,t)-\delta_{C}(z,z)\\}\widetilde{\psi}(z,t)\right]$
(8)

and

$\displaystyle\frac{\partial}{\partial t}\widetilde{\psi}^{+}(z,t)$
$\displaystyle=+\frac{i}{\hslash}\left[-\frac{\hslash^{2}}{2m}\frac{\partial^{2}}{\partial
z^{2}}\widetilde{\psi}^{+}(z,t)+V(z,t)\widetilde{\psi}^{+}(z,t)+g\\{\widetilde{\psi}^{+}(z,t)\widetilde{\psi}(z,t)-.\delta_{C}(z,z)\\}\widetilde{\psi}^{+}(z,t)\right]$
(9)

where the $\widetilde{\psi}(z,t),\widetilde{\psi}^{+}(z,t)$ are time-dependent
stochastic field functions that now replace the original phase-space field
functions $\psi(z),\psi^{+}(z)$. These are functions of position which may be
expanded in terms of an independent set of (non-stochastic) basis functions or
modes and where the amplitudes (or stochastic phase variables) for the mode
expansion are treated as stochastic quantities. Stochasticity in the context
of phase space theory is discussed in Appendix G of Ref. [22]. As there are no
diffusion terms in the FFPE, there are no Gaussian-Markov noise terms
satisfying the standard results $\overline{\Gamma_{a}\,(t)}=0$,
$\overline{\Gamma_{a}\,(t_{1})\Gamma_{b}\,(t_{2})}$
$=\delta_{a,b}\delta(t_{1}-t_{2})$ etc, in the Ito SFE for the truncated
Wigner $W^{+}$ treatment (the bar denotes a stochastic average). The solutions
to the Ito SFE are uniquely determined from the initial functions
$\widetilde{\psi}(z,0)$ and $\widetilde{\psi}^{+}(z,0)$, as if the fields were
like classical fields. However,
$\widetilde{\psi}(z,t),\widetilde{\psi}^{+}(z,t)$ are non-classical and
stochastic because the initial functions $\widetilde{\psi}(z,0)$ and
$\widetilde{\psi}^{+}(z,0)$ are members of an ensemble of stochastic fields,
with a distribution chosen to represent the features of the initial quantum
state $\widehat{\rho}(0)$. This in turn leads to an ensemble of time-dependent
stochastic fields $\widetilde{\psi}(z,t),\widetilde{\psi}^{+}(z,t)$
representing the time evolution of the quantum state, this being equivalent to
the time evolution in the Wigner distribution functional
$W[\mathbf{\psi}(z)].$

####  2.2.3 Mean-Field Theory and Gross-Pitaevski Equation

The Ito stochastic field equations resemble the _Gross-Pitaevskii equation_
(GPE) that applies for the mean-field approximation in which all bosons are
assumed to occupy a single condensate mode. The GPE would represent another
approach to treating the periodically driven BEC, but as it does not allow for
the presence of other modes this mean-field approach cannot allow for quantum
fluctuation effects that are treated via the TWA. The time-dependent GPE for
the condensate wave function $\Phi_{c}(z,t)$ \- normalised as
${\displaystyle\int}dz\,\Phi_{c}(z,t)^{\ast}$ $\Phi_{c}(z,t)$ $=N$ \- is given
by

$\frac{\partial}{\partial
t}\Phi_{c}(z,t)=-\frac{i}{\hslash}\left[-\frac{\hslash^{2}}{2m}\frac{\partial^{2}}{\partial
z^{2}}\Phi_{c}(z,t)+V(z,t)\Phi_{c}(z,t)+g\\{\Phi_{c}(z,t)^{\ast}\Phi_{c}(z,t)\\}\Phi_{c}(z,t)\right]$
(10)

If we choose $\widetilde{\psi}^{+}(z,t)=\widetilde{\psi}(z,t)^{\ast}$ and
$\widetilde{\psi}(z,t)=\Phi_{c}(z,t)$ the Ito SFE has almost the same form as
the GPE. However, as well as $\Phi_{c}(z,t)$ being non-stochastic, the
quantity $\delta_{C}(z,z)$ is absent from the non-linear term in the GPE that
allows for boson-boson interactions in the mean-field approximation. As
${\displaystyle\int}dz\,\Phi_{c}(z,t)^{\ast}\Phi_{c}(z,t)=N$ and
${\displaystyle\int}dz\,\delta_{C}(z,z)=n_{M}$ (where $n_{M}$ is the number of
modes that need to be used to treat the physics), we see that the
$\delta_{C}(z,z)$ term in the Ito SFE will be relatively unimportant compared
to the $\widetilde{\psi}^{+}(z,t)\widetilde{\psi}(z,t)$ term in the usual case
for a BEC system, where the number of bosons considerably exceeds the number
of modes that need to be considered. The stochasticity of the initial
$\widetilde{\psi}(z,0)$ and $\widetilde{\psi}^{+}(z,0)$ should be more
important than the $\delta_{C}(z,z)$ term in treating quantum effects for
interacting bosons.

We define throughout a unity normalised _condensate mode function_
$\psi_{c}(z,t)$ as

$\psi_{c}(z,t)=\Phi_{c}(z,t)/\sqrt{N}$ (11)

####  2.2.4 Validity and Reliability of Truncated Wigner Approximation

Since mean-field theory or time-dependent Bogoliubov theory can be derived
from the TWA as approximations (as shown in Appendix 6), we expect the TWA to
be more reliable than both of these approximate approaches - which themselves
have been successfully used to treat cold bosonic gases. Indeed, the TWA has
already been extensively used to describe cold bosonic gases. However, as
explained above, the TWA approach involves a key approximation – the neglect
of third-order derivative terms in the FFPE. Ultimately, the reliability of
the TWA rests on whether these terms can be neglected, but as far as we are
aware no calculations of the size of these terms has been carried out. So, at
present, the only tests of the reliabilty of the TWA are: (a) whether its
predictions agree with experiment - and the relevant experiment has not yet
been carried out; and/or (b) whether it agrees with exact full FFPE
calculations – and these have not been done. Detailed issues are discussed in
Appendix 7.

####  2.2.5 Mean Energy

The expression for the _mean energy_

$\left\langle\widehat{H}\right\rangle=Tr(\widehat{H}\,\widehat{\rho})$ (12)

in terms of the stochastic field functions
$\widetilde{\psi}(z,t),\widetilde{\psi}^{+}(z,t)$ is

$\displaystyle\left\langle\widehat{H}\right\rangle$
$\displaystyle={\displaystyle\int}dz\,\left(\frac{\hbar^{2}}{2m}\overline{\frac{\partial}{\partial
z}\widetilde{\psi}^{+}(z,t)\,\frac{\partial}{\partial
z}\widetilde{\psi}(z,t)}\right)\,$
$\displaystyle+{\displaystyle\int}dz\,\left(V(z,t)\,\overline{\widetilde{\psi}^{+}(z,t)\,\,\widetilde{\psi}(z,t)}\right)$
$\displaystyle+{\displaystyle\int}dz\,\left(\frac{g}{2}\;\overline{\widetilde{\psi}^{+}(z,t)\,\widetilde{\psi}^{+}(z,t)\,\widetilde{\psi}(z,t)\,\widetilde{\psi}(z,t)}\right)$
$\displaystyle-{\displaystyle\int}dz\,\left(g\,\delta_{C}(z,z)\;\overline{\widetilde{\psi}^{+}(z,t)\,\,\widetilde{\psi}(z,t)}\right)$
$\displaystyle-{\displaystyle\int}dz\,\left(\frac{\hbar^{2}}{4m}(\triangle\delta_{C}(z,z))\right)-{\displaystyle\int}dz\,\left(\frac{1}{2}\delta_{C}(z,z)\;V(z,t)\right)$
$\displaystyle+{\displaystyle\int}dz\,\left(\frac{g}{4}\,\delta_{C}(z,z)^{2}\right)$
(13)

where

$\triangle\delta_{C}(z,z)={\displaystyle\sum\limits_{k}}\left(\frac{\partial}{\partial
z}\phi_{k}(z)\right)\left(\frac{\partial}{\partial
z}\phi_{k}(z)\right)^{\ast}$ (14)

and $\phi_{k}(z)$ being any set of suitable orthonormal mode functions. The
proof is given in Appendix 8.

For the present situation where the number of bosons considerably exceeds the
number of modes that need to be considered, only the first three terms are
important. The last three terms are non-stochastic. An expression for the mean
energy is given below in Eq. (51) in terms of gravitational modes.

####  2.2.6 Position Probability Density

The _number operator_ is given by

$\widehat{N}={\displaystyle\int}dz\,\hat{\Psi}^{{\dagger}}(z)\hat{\Psi}(z)$
(15)

which is the integral over an operator
$\widehat{N}(z)=\hat{\Psi}^{{\dagger}}(z)\hat{\Psi}(z)$ associated with the
position probability density for the bosons in the system. We consider states
which are eigenstates of $\widehat{N}$ with eigenvalue $N$.

The _position probability density_ (PPD) is defined as the mean value of
$\widehat{N}(z)$

$F(z,t)=Tr(\hat{\Psi}^{{\dagger}}(z)\hat{\Psi}(z)\rho(t))$ (16)

This essentially specifies the relative numbers of bosons found at various
positions $z$ as a function of $t$.

As indicated in the Introduction, to employ the $W^{+}$ distribution we must
express $\hat{\Psi}^{{\dagger}}(z)\hat{\Psi}(z)$ in the symmetric form
$\\{\hat{\Psi}^{{\dagger}}(z)\hat{\Psi}(z)\\}=\frac{1}{2}(\hat{\Psi}^{{\dagger}}(z)\hat{\Psi}(z)+\hat{\Psi}(z)\hat{\Psi}^{{\dagger}}(z))$
via
$\hat{\Psi}^{{\dagger}}(z)\hat{\Psi}(z)=\\{\hat{\Psi}^{{\dagger}}(z)\hat{\Psi}(z)\\}-\frac{1}{2}\delta_{C}(z,z)$,
so that the position number density is then given either in terms of a
functional integral or equivalently via a stochastic average

$\displaystyle F(z,t)$
$\displaystyle={\textstyle\int}{\textstyle\int}D^{2}\psi
D^{2}\psi^{+}\,\psi(z)W[\mathbf{\psi}(z)]\psi^{+}(z)-\frac{1}{2}\delta_{C}(z,z)$
(17)
$\displaystyle=\overline{\widetilde{\psi}(z,t)\widetilde{\psi}^{+}(z,t)}-\frac{1}{2}\delta_{C}(z,z)$
(18)

The stochastic average is the expression used in numerical calculations. Note
the presence of the $\frac{1}{2}\delta_{C}(z,z)$, but as the spatial integral
of this term is of order $\frac{1}{2}n_{M}$, whilst that of
$\overline{\widetilde{\psi}(z,t)\widetilde{\psi}^{+}(z,t)}$ is of order
$N_{c}\gg n_{M}$, the second term in Eq. (18) should not be important.

It should be noted that other distribution functionals can also be used to
represent the density operator, and these have their own functional Fokker-
Planck equations and Ito stochastic field equations. These include the
positive $P+$ distribution functional, where here the FFPE includes a
diffusion term and the Ito SFE involve Gaussian-Markoff noise terms. The $P+$
case is discussed in Ref [22] (see Sections 15.1.5 and 15.1.7).

A direct numerical solution of the Ito stochastic field equations for each
initial $\widetilde{\psi}(z,0)$ and $\widetilde{\psi}^{+}(z,0)$ is one
possible approach to calculating the position probability density via Eq (18),
with an ensemble of different initial stochastic fields chosen to represent
the properties of the initial quantum state. In such a treatment the
underlying presence of mode functions and their frequencies is not made
explicit. Our approach however will involve mode expansions.

####  2.2.7 Quantum Correlation Function

As well as the position probability density, there is a further quantity that
is of interest in describing BECs. This is the first-order _quantum
correlation function_ (QCF) which is defined [28] as

$P(z,z^{\\#},t)=Tr(\hat{\Psi}^{{\dagger}}(z^{\\#})\hat{\Psi}(z)\rho(t))$ (19)

where $z,z^{\\#}$ are two spatial positions. Since
$\hat{\Psi}^{{\dagger}}(z^{\\#})\hat{\Psi}(z)=\\{\hat{\Psi}^{{\dagger}}(z^{\\#})\hat{\Psi}(z)\\}-\frac{1}{2}\delta_{C}(z,z^{\\#})$,
the QCF can also be expressed as a stochastic average via

$P(z,z^{\\#},t)=\overline{\widetilde{\psi}(z,t)\widetilde{\psi}^{+}(z^{\\#},t)}-\frac{1}{2}\delta_{C}(z,z^{\\#})$
(20)

This QCF can be expressed in terms of _natural orbitals_ $\chi_{i}(z,t)$ and
their _occupation numbers_ $p_{i}$, as will now be shown. First, we see that
$P(z,z^{\\#},t)=P(z^{\\#},z,t)^{\ast}$. Then if we introduce an orthonormal
set of mode functions $\phi_{k}(z)$ we can define a matrix $P$ with elements
$P_{k,l}={\displaystyle\iint}dz_{1}\,dz_{1}^{\\#}\,\phi_{k}^{\ast}(z)\,P(z,z^{\\#},t)\,\phi_{l}(z^{\\#})$.
It is easy to show that $P$ is Hermitian, so it can be diagonalised via a
unitary matrix $U$ such that $P=U\Delta U^{{\dagger}}$, where $\Delta$ is a
diagonal matrix containing the eigenvalues $p_{i}$ of $P$. These eigenvalues
are real and positive. Substituting for $P$ we find that the QCF can be
written in terms of natural orbitals $\chi_{i}(z,t)$ and their occupation
numbers $p_{i}$ (see Eq.(21)) as

$P(z,z^{\\#},t)={\textstyle\sum\limits_{i=0}}p_{i}\,\chi_{i}(z,t)\,\chi_{i}(z^{\\#},t)^{\ast}$
(21)

where the natural orbitals are

$\chi_{i}(z,t)={\textstyle\sum\limits_{k}}U_{k,i}\,\phi_{k}(z,t)$ (22)

Situations where almost all bosons occupy one natural orbital are an indicator
of a BEC. By convention the natural orbital with the largest occupancy is
listed as $i=0$. The QCF essentially specifies the long-range spatial
coherence that applies in BECs.

####  2.2.8 One-Body Projector

We can define a _one-body projection_ operator (OBP) onto any single particle
state $\left|\phi\right\rangle$ for an $n$ body system in first quantisation
as

$\widehat{M}={\displaystyle\sum\limits_{i=1}^{n}}(\left|\phi\right\rangle\left\langle\phi\right|)_{i}$
(23)

the sum being over _all_ the identical bosons. Note that for a system of
identical bosons the operator must be _symmetric_ under particle permutations.
Typically in a BEC with $N\gg 1$ bosons in the condensate mode the single-
particle state would be taken as the initial condensate mode function with
$\left\langle z|\psi_{c}\right\rangle=$ $\psi_{c}(z,0)$. This quantity is a
Hermitian operator for the many-body system and can be regarded as an
observable. Its mean value could therefore be measured as a function of time.

The operator $\widehat{M}$ can be written in second quantisation form using an
orthonormal set of mode functions $\phi_{k}(z)$ and their annihilation,
creation operators $\widehat{a}_{k}$, $\widehat{a}_{k}^{{\dagger}}$ using
standard procedures in which the field operators are expanded in terms of a
set of orthonormal mode functions $\phi_{k}(z)$ as
$\widehat{\Psi}(z)={\displaystyle\sum\limits_{k}}\phi_{k}(z)\,\widehat{a}_{k}$
and
$\widehat{\Psi}(z)^{{\dagger}}={\displaystyle\sum\limits_{k}}\phi_{k}(z)^{\ast}\,\widehat{a}_{k}^{{\dagger}}$.
We have for the case of the condensate mode function

$\displaystyle\widehat{M}_{c}$
$\displaystyle={\displaystyle\sum\limits_{k,l}}\widehat{a}_{k}^{{\dagger}}\,\widehat{a}_{l}\,\left\langle\phi_{k}(1)|\,(\left|\psi_{c}\right\rangle\left\langle\psi_{c}\right|)_{1}\,|\,\phi_{l}(1)\right\rangle$
$\displaystyle={\displaystyle\sum\limits_{k,l}}\widehat{a}_{k}^{{\dagger}}\,\widehat{a}_{l}\,{\displaystyle\int}dz_{1}\,\,\phi_{k}^{\ast}(z_{1})\psi_{c}(z_{1},0){\displaystyle\int}dz_{1}^{\\#}\,\,\psi_{c}^{\ast}(z_{1}^{\\#},0)\phi_{l}(z_{1}^{\\#})$
$\displaystyle={\displaystyle\sum\limits_{k,l}}\left({\displaystyle\int}dz^{\\#}\,\phi_{k}(z^{\\#})\widehat{\Psi}(z^{\\#})^{{\dagger}}\right)\left({\displaystyle\int}dz\,\phi_{l}^{\ast}(z)\widehat{\Psi}(z)\right)\left({\displaystyle\int}{\displaystyle\int}dz_{1}\,dz_{1}^{\\#}\,\phi_{k}^{\ast}(z_{1})\psi_{c}(z_{1},0)\psi_{c}^{\ast}(z_{1}^{\\#},0)\phi_{l}(z_{1}^{\\#})\right)$
$\displaystyle={\displaystyle\int}{\displaystyle\int}dz\,dz^{\\#}\,\psi_{c}(z^{\\#},0)\,\widehat{\Psi}(z^{\\#})^{{\dagger}}\,\widehat{\Psi}(z)\,\psi_{c}^{\ast}(z,0)$
(24)

where the completeness results
${\displaystyle\sum\limits_{k}}\phi_{k}(z^{\\#})\,\phi_{k}^{\ast}(z_{1})=\delta(z^{\\#}-z_{1})$
and
${\displaystyle\sum\limits_{k}}\phi_{l}(z_{1}^{\\#})\,\phi_{l}^{\ast}(z)=\delta(z-z_{1}^{\\#})$
are used and the integrals over $z_{1}$ and $z_{1}^{\\#}$ are carried out.
Thus, the one-body projector operator involves an integral of the position
correlation operator $\widehat{\Psi}(z^{\\#})^{{\dagger}}\,\widehat{\Psi}(z)$
times $\psi_{c}(z^{\\#},0)\psi_{c}^{\ast}(z,0)$.

Hence, we see that the mean value for the one-body projector operator is

$M_{c}(t)={\displaystyle\int}{\displaystyle\int}dz\,dz^{\\#}\,\psi_{c}(z^{\\#},0)\,\,\psi_{c}^{\ast}(z,0)\,P(z,z^{\\#},t)$
(25)

which is the double space integral of
$\psi_{c}(z^{\\#},0)\psi_{c}^{\ast}(z,0)$ multiplied by the time dependent
_first order quantum correlation function_ $P(z,z^{\\#},t)$, which is defined
and its expression given in terms of stochastic field functions as in Eqs.
(19) and (20).

If we substitute for the QCF in terms of the stochastic field functions
$\widetilde{\psi}(z,t)$ and $\widetilde{\psi}^{+}(z,t)$ we get

$M_{c}(t)={\displaystyle\int}{\displaystyle\int}dz\,dz^{\\#}\,\psi_{c}(z^{\\#},0)\,\psi_{c}^{\ast}(z,0)\left\\{\overline{\widetilde{\psi}(z,t)\widetilde{\psi}^{+}(z^{\\#},t)}-\frac{1}{2}\delta_{C}(z,z^{\\#})\right\\}$
(26)

Note that the result involves the condensate mode function $\psi_{c}(z,0)$ at
time $0$, rather than at time $t$, as in the expression in Eq. (39) for the
occupancy $N_{C}(t)$ of the time-dependent condensate mode. Thus, the mean
value of the one-body projector is a different quantity to the quantum
depletion - defined below in Section 2.2.10. We can use the mean value of the
one-body projector $M_{c}(t)$ as an observable. This can obviously be
calculated from the time-dependent stochastic field functions
$\widetilde{\psi}(z,t)$ and $\widetilde{\psi}^{+}(z^{\\#},t)$ and the
condensate mode function $\psi_{c}(z,0)$ as a function of time. The term
$\frac{1}{2}\delta_{C}(z,z^{\\#})$ would be small compared to the
$\overline{\widetilde{\psi}(z,t)\widetilde{\psi}^{+}(z^{\\#},t)}$ term for
$N_{c}\gg 1$. The quantity $M_{c}(t)$ has the desired feature of revealing the
periodicity of the QCF via its Fourier transform (FT).

Also we have

$M_{c}(0)=N$ (27)

assuming the initial condition of all bosons in mode $\psi_{c}$ and using Eq.
(89) (see below). It is straightforward to rephrase the above treatment in a
normalised form $M_{c}(t)/M_{c}(0)$ using Eq. (27) for the initial value of
$M_{c}(t)$

The Fourier transform of the position probability density reveals the
frequencies at which the PPD is oscillating. By taking the FT of $M_{c}(t)$ we
have

$FT\\{M_{c}(t)\\}={\displaystyle\int}{\displaystyle\int}dz\,dz^{\\#}\,\psi_{c}(z^{\\#},0)\,\,\psi_{c}^{\ast}(z,0)\,\;FT\\{P(z,z^{\\#},t)\\}$
(28)

since the initial condensate mode function $\,\psi_{c}(z,0)$ is time-
independent. Thus the FT of $M_{c}(t)$ would reveal the _same frequencies_ as
the FT of the QCF.

Hence, the one-body projector approach yields a useful indicator for
describing the behaviour of the QCF (and also of its diagonal terms when
$z=z^{\\#}$ \- the position probability density), whose periodicities
determine whether time crystal behaviour is occurring. The FT of $M_{c}(t)$
does not of course explain the periodicity of the QCF or the PPD - that
explanation requires a consideration of the numerous time scales involved in
the system (the drive period $T$, the gravitational mode frequencies, the time
scale associated with boson-boson interaction, the bounce time of the BEC,
etc.).

The mean value of the one body projector $M_{c}(t)$ for a many-body system is
similar to the autocorrelation function or the _fidelity_ $F_{K}(t)$
introduced by Kuros et. al. [19] (see Eq. (6) therein) in a mean-field theory.
The fidelity was defined by

$F_{K}(t)=|\left\langle\Phi_{c}(0)|\Phi_{c}(t)\right\rangle|^{2}=|{\displaystyle\int}dz\;\Phi_{c}(z,0)\Phi_{c}^{\ast}(z,t)\,|^{2}$
(29)

where $\Phi_{c}(z,t)$ is the condensate wave-function, given by the solution
of the GPE. This specifies how similar the time-dependent condensate wave-
function is to its initial form.

This can be written as (note the position orders)

$F_{K}(t)={\displaystyle\int}{\displaystyle\int}dz\,dz^{\\#}\,P_{K}(z^{\\#},z,0)\,P_{K}(z,z^{\\#},t)$
(30)

where

$P_{K}(z,z^{\\#},t)\,=\Phi_{c}(z,t)\Phi_{c}^{\ast}(z^{\\#},t)$ (31)

is the QCF in the single-mode mean-field approximation. The similarity of
Eq.(30) to Eq.(25) for the mean value of the condensate mode function OBP
$M_{c}(t)$ is clear. The Fourier transform of $F_{K}(t)$ is related to the FT
of the QCF $P_{K}(z,z^{\\#},t)$ via a similar relationship to that in Eq.
(28). Note that $F_{K}(t)$ is not proportional to
${\displaystyle\int}dz\,\rho(z,0)\,\rho(z,t)$ \- as stated in Ref. [19], where
$\rho(z,t)=|\Phi_{c}(z,t)|^{2}$ is the mean-field position probability
density.

####  2.2.9 Other One-Body Projectors

The one-body projector concept can be extended to cases where the mode
function chosen is different to the condensate mode $\psi_{c}(z,0).$ One case
of particular interest is where the mode function is one of the two Wannier-
like modes $\Phi_{i}(z,t)$ $(i=1,2)$ (see Sect 2.3.3 and Eq (59) below).

In this case we would have for the OBP operator $\widehat{M}_{Wi}(t)$ and its
mean value $M_{Wi}(t)\equiv N_{i}(t)$

$\displaystyle\widehat{M}_{Wi}(t)$
$\displaystyle={\displaystyle\int}{\displaystyle\int}dz\,dz^{\\#}\,\Phi_{i}(z^{\\#},t)\,\widehat{\Psi}(z^{\\#})^{{\dagger}}\,\widehat{\Psi}(z)\,\Phi_{i}^{\ast}(z,t)$
$\displaystyle M_{Wi}(t)$
$\displaystyle=Tr\left(\widehat{M}_{Wi}(t)\widehat{\rho}(t)\right)\equiv
N_{i}(t)$ $\displaystyle N_{i}(t)$
$\displaystyle={\displaystyle\int}{\displaystyle\int}dz\,dz^{\\#}\,\Phi_{i}(z^{\\#},t)\,\,\Phi_{i}^{\ast}(z,t)\left\\{\overline{\widetilde{\psi}(z,t)\widetilde{\psi}^{+}(z^{\\#},t)}-\frac{1}{2}\delta_{C}(z,z^{\\#})\right\\}\qquad
i=1,2$ (32)

####  2.2.10 Quantum Depletion from Condensate Mode

We can expand the field operators, field functions and the stochastic field
functions in terms of a time-dependent condensate mode function
$\psi_{c}(z,t)=\Phi_{c}(z,t)/\sqrt{N}$ obtained from the solution of the time-
dependent GPE (10) and any set of orthogonal non-condensate modes
$\psi_{k}(z,t)$ as

$\displaystyle\widehat{\Psi}(z)$
$\displaystyle=\psi_{c}(z,t)\widehat{c}_{0}(t)+{\displaystyle\sum\limits_{k\neq
0}}\psi_{k}(z,t)\widehat{c}_{k}(t)\qquad\widehat{\Psi}(z)^{{\dagger}}=\psi_{c}(z,t)^{\ast}\widehat{c}_{0}^{{\dagger}}(t)+{\displaystyle\sum\limits_{k\neq
0}}\psi_{k}(z,t)^{\ast}\widehat{c}_{k}^{{\dagger}}(t)$ (33)
$\displaystyle\widetilde{\psi}(z,t)$
$\displaystyle=\psi_{c}(z,t)\widetilde{\gamma}_{0}(t)+{\displaystyle\sum\limits_{k\neq
0}}\psi_{k}(z,t)\widetilde{\gamma}_{k}(t)\qquad\widetilde{\psi}(z,t)^{+}=\psi_{c}(z,t)^{\ast}\widetilde{\gamma}_{0}^{+}(t)+{\displaystyle\sum\limits_{k\neq
0}}\psi_{k}(z,t)^{\ast}\widetilde{\gamma}_{k}^{+}(t)$ (34)

where $\widehat{c}_{0}(t)$, $\widehat{c}_{0}^{{\dagger}}(t)$ are condensate
mode annihilation, creation operators and $\widetilde{\gamma}_{0}(t)$,
$\widetilde{\gamma}_{0}^{+}(t)$ are the related stochastic amplitudes. Similar
operators $\widehat{c}_{k}(t)$, $\widehat{c}_{k}^{{\dagger}}(t)$ and
stochastic amplitudes $\widetilde{\gamma}_{k}(t)$,
$\widetilde{\gamma}_{k}^{+}(t)$ apply for the non-condensate modes.

The number of bosons left in the original condensate mode $N_{C}(t)$ is given
by

$N_{C}(t)=\left\langle\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}\right\rangle.$
(35)

We can now derive an expression for the quantum depletion (QD)

$N_{D}(t)=N-N_{C}(t)$ (36)

which specifies the loss of bosons from the condensate mode and is another
physical quantity of interest. The quantum depletion gives a measure of the
failure of the mean-field theory

As
$\\{\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}\\}=\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}+1/2$
we can use the Wigner distribution approach to replace
$\left\langle\\{\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}\\}\right\rangle$ by
$\overline{\widetilde{\gamma}_{0}^{+}(t)\widetilde{\gamma}_{0}(t)}$ so that we
have

$N_{C}(t)=\overline{\widetilde{\gamma}_{0}^{+}(t)\widetilde{\gamma}_{0}(t)}-1/2$
(37)

By using Eq. (83) for the stochastic fields $\widetilde{\Psi}(z,t)$,
$\widetilde{\Psi}^{+}(z,t)$ the stochastic phase amplitudes for the condensate
mode are given by

$\widetilde{\gamma}_{0}(t)=\int
dz\,\psi_{c}^{\ast}(z,t)\widetilde{\Psi}(z,t)\qquad\widetilde{\gamma}_{0}^{+}(t)=\int
dz^{\\#}\,\psi_{c}(z^{\\#},t)\widetilde{\Psi}^{+}(z^{\\#},t)$ (38)

Hence, using Eqs. (19) and (20) and
$\delta_{C}(z,z^{\\#})=\psi_{c}(z,t)\psi_{c}(z^{\\#},t)^{\ast}+{\displaystyle\sum\limits_{k\neq
c}}\psi_{k}(z,t)\psi_{k}(z^{\\#},t)^{\ast}$ we can then show that

$\displaystyle N_{C}(t)$
$\displaystyle={\displaystyle\int}{\displaystyle\int}dz\,dz^{\\#}\,\psi_{c}(z^{\\#},t)\left\langle\widehat{\Psi}(z^{\\#})^{{\dagger}}\,\widehat{\Psi}(z)\right\rangle\psi_{c}(z,t)^{\ast}$
$\displaystyle={\displaystyle\int}{\displaystyle\int}dz\,dz^{\\#}\,\psi_{c}(z^{\\#},t)\psi_{c}(z,t)^{\ast}P(z,z^{\\#},t)$
(39)

where $P(z,z^{\\#},t)$ is the first-order quantum correlation function. Noting
the result (25) for the mean value of the one-body projector, we see that the
quantum depletion has a similar relation to the QCF, but now involving the
time-dependent condensate mode function. Taking the Fourier transform of each
side of Eq. (39) results in the FT of the quantum depletion not reflecting the
same periodicities as in the QCF due to the presence of the
$\psi_{c}(z^{\\#},t)\psi_{c}(z,t)^{\ast}$factor. However, as noted above the
quantum depletion has a different role, namely indicating whether or not the
mean-field theory still applies.

By expressing $P(z,z^{\\#},t)$ in terms of natural orbitals via (21) we see
that the condensate mode occupancy is also given by

$N_{C}(t)={\displaystyle\sum\limits_{i=0}}p_{i}\times\left|{\displaystyle\int}dz\,\psi_{c}^{\ast}(z,t)\chi_{i}(z,t)\right|^{2}$
(40)

Subsituting for the QCF in terms of stochastic fields gives

$\displaystyle N_{C}(t)$
$\displaystyle={\displaystyle\int}{\displaystyle\int}dz\,dz^{\\#}\,\psi_{c}(z^{\\#},t)\psi_{c}(z,t)^{\ast}\left(\overline{\widetilde{\psi}(z,t)\widetilde{\psi}^{+}(z^{\\#},t)}-\frac{1}{2}\delta_{C}(z,z^{\\#})\right)$
$\displaystyle={\displaystyle\int}{\displaystyle\int}dz\,dz^{\\#}\,\psi_{c}(z^{\\#},t)\psi_{c}(z,t)^{\ast}\left(\overline{\widetilde{\psi}(z,t)\widetilde{\psi}^{+}(z^{\\#},t)}\right)\;-\frac{1}{2}$
(41)

where the result
$\delta_{C}(z,z^{\\#})=\psi_{c}(z,t)\psi_{c}(z^{\\#},t)^{\ast}+{\displaystyle\sum\limits_{k\neq
c}}\psi_{k}(z,t)\psi_{k}(z^{\\#},t)^{\ast}$ has again been used.

###  2.3 Periodic Driven BEC- Stochastic Phase Amplitudes for Separate Modes

In this section we set out the dynamical equations for the stochastic phase-
space amplitudes associated with separate modes that are involved when the
field operator, field function and stochastic field are expanded in terms of a
suitable set of orthogonal mode functions. These stochastic amplitude
equations describe the time evolution of the periodically driven BEC. Later,
we will consider other mode functions (Bogoliubov modes) that are more
suitable for describing the preparation of the BEC in a trap potential before
it is released and subjected to periodic driving. There are, however, several
alternative choices that are all suitable for treating the driven BEC.

####  2.3.1 Floquet Modes

As we are considering the effect of a periodic driving field on the BEC one
possibility is to choose Floquet modes [29], [14], which are essentially
single boson wave functions for evolution due to the periodic field. They
satisfy the equation

$-\frac{\hslash^{2}}{2m}\frac{\partial^{2}}{\partial
z^{2}}\phi_{k}(z,t)+V(z,t)\phi_{k}(z,t)-i\hbar\frac{\partial}{\partial
t}\phi_{k}(z,t)=\hslash\nu_{k}\phi_{k}(z,t)$ (42)

and are periodic $\phi_{k}(z,t)=\phi_{k}(z,t+T)$ as well as being orthonormal
${\displaystyle\int}dz\,\phi_{k}(z,t)^{\ast}\phi_{l}(z,t)=\delta_{k,l}$ at all
times $t$. The Floquet frequencies $\nu_{k}$ are time independent, and form
zones analogous to the Brillouin zones that occur for particles moving in a
spatially periodic potential.

Thus, we have

$\displaystyle\widehat{\psi}(z)$
$\displaystyle={\displaystyle\sum\limits_{k}}\widehat{a}_{k}(t)\,\phi_{k}(z,t)\qquad\widehat{\psi}^{{\dagger}}(z)={\displaystyle\sum\limits_{k}}\widehat{a}_{k}^{{\dagger}}(t)\,\phi_{k}^{\ast}(z,t)$
(43) $\displaystyle\psi(z)$
$\displaystyle={\displaystyle\sum\limits_{k}}\alpha_{k}(t)\,\phi_{k}(z,t)\qquad\psi^{+}(z)={\displaystyle\sum\limits_{k}}\alpha_{k}^{+}(t)\,\phi_{k}^{\ast}(z,t)$
(44) $\displaystyle\widetilde{\psi}(z,t)$
$\displaystyle={\displaystyle\sum\limits_{k}}\widetilde{\alpha}_{k}(t)\,\phi_{k}(z,t)\qquad\widetilde{\psi}^{+}(z,t)={\displaystyle\sum\limits_{k}}\widetilde{\alpha}_{k}^{+}(t)\,\phi_{k}^{\ast}(z,t)$
(45)

are the expansions for the field operators, field functions and stochastic
field that occur in the Hamiltonian, the FFPE and the Ito SFE, respectively.
The mode annihilation, creation operators for the Floquet modes are time
dependent, but still satisfy the standard Bose commutation rules
$[\widehat{a}_{k},\widehat{a}_{l}]=$
$[\widehat{a}_{k}^{{\dagger}},\widehat{a}_{l}^{{\dagger}}]=0,$
$[\widehat{a}_{k},\widehat{a}_{l}^{{\dagger}}]=\delta_{k,l}$. The phase-space
variables $\alpha_{k},\alpha_{k}^{+}$ that represent the mode annihilation,
creation operators are also time dependent - as are their related stochastic
phase-space variables $\widetilde{\alpha}_{k},\widetilde{\alpha}_{k}^{+}$.

However, using Floquet modes in the numerics has the disadvantage that they
are time dependent, so sets of such modes must be stored on the computer as
functions of time. Instead the numerics are more conveniently carried out
using time-independent gravitational modes. For completeness, we have set out
the equations that would be used for the evolution of stochastic phase
variables for Floquet modes together with Floquet-based expressions for the
position probability density and QCF in Appendix 9.

####  2.3.2 Gravitational Modes

Rather than Floquet modes - since we are considering the effect of a static
gravitational field on the BEC - it is convenient to choose gravitational
modes, which are essentially single boson wave-functions for evolution due to
the static gravitational field. These have the advantage of being time
independent as well as being orthogonal. They satisfy the equation

$-\frac{\hslash^{2}}{2m}\frac{\partial^{2}}{\partial
z^{2}}\xi_{k}(z)+mg_{E}z\,\xi_{k}(z)=\hslash\epsilon_{k}\xi_{k}(z)$ (46)

Expansion of the field function and stochastic field function in terms of
these modes gives

$\displaystyle\psi(z)$
$\displaystyle={\displaystyle\sum\limits_{k}}\eta_{k}\,\xi_{k}(z)\qquad\psi^{+}(z)={\displaystyle\sum\limits_{k}}\eta_{k}^{+}\,\xi_{k}^{\ast}(z)$
(47) $\displaystyle\widetilde{\psi}(z,t)$
$\displaystyle={\displaystyle\sum\limits_{k}}\widetilde{\eta}_{k}(t)\,\xi_{k}(z)\qquad\widetilde{\psi}^{+}(z,t)={\displaystyle\sum\limits_{k}}\widetilde{\eta}_{k}^{+}(t)\,\xi_{k}^{\ast}(z)$
(48)

The phase-space variables $\eta_{k},\eta_{k}^{+}$ that represent the mode
annihilation, creation operators are time independent - whereas their related
stochastic phase-space variables
$\widetilde{\eta}_{k},\widetilde{\eta}_{k}^{+}$ are time dependent.

##### Position ProbabilityDensity, QCF and Mean Energy - Gravitational Modes

The position probability density in Eq. (18) can be expressed in terms of
gravitational mode functions as

$F(z,t)={\displaystyle\sum\limits_{k,l}}\xi_{k}(z)\,\xi_{l}^{\ast}(z)\,\left[\overline{\widetilde{\eta}_{k}(t)\,\widetilde{\eta}_{l}^{+}(t)}\,-\frac{1}{2}\delta_{k,l}\right]$
(49)

and involves the stochastic average of products of stochastic phase-space
variables.

Similarly, the QCF in Eq. (20) can also be expressed in terms of gravitational
mode functions (see Eq.(21)) as

$P(z,z^{\\#},t)={\displaystyle\sum\limits_{k,l}}\xi_{k}(z)\,\xi_{l}^{\ast}(z^{\\#})\,\left[\overline{\widetilde{\eta}_{k}(t)\,\widetilde{\eta}_{l}^{+}(t)}\,-\frac{1}{2}\delta_{k,l}\right]$
(50)

Furthermore, the mean energy in Eq. (13) can be simplified by using the
expansion of the stochastic field functions (48) in terms of gravitational
modes along with the defining Eq. (46) for these modes. We have

$\displaystyle\left\langle\widehat{H}\right\rangle$
$\displaystyle={\displaystyle\sum\limits_{k=0}}\hslash\epsilon_{k}\,\left(\overline{\widetilde{\eta}_{k}(t)\,\widetilde{\eta}_{k}^{+}(t)}\,-\frac{1}{2}\right)$
$\displaystyle-mg_{E}\lambda\cos\omega
t\,{\displaystyle\int}dz\,z\,\left(\overline{\widetilde{\psi}^{+}(z,t)\,\widetilde{\psi}(z,t)}-\frac{1}{2}\delta_{C}(z,z)\right)$
$\displaystyle+\frac{g}{2}{\displaystyle\int}dz\,\left(\overline{\widetilde{\psi}^{+}(z,t)\,\widetilde{\psi}^{+}(z,t)\,\widetilde{\psi}(z,t)\,\widetilde{\psi}(z,t)}-2\delta_{C}(z,z)\,\overline{\widetilde{\psi}^{+}(z,t)\,\widetilde{\psi}(z,t)}+\frac{1}{2}\delta_{C}(z,z)^{2}\right)$
(51)

The derivation is set out in Appendix 8.

##### Evolution of Stochastic Phase Variables for Gravitational Modes

Coupled equations for the stochastic phase-space variables
$\widetilde{\eta}_{k},\widetilde{\eta}_{k}^{+}$ can be obtained allowing for
the effect of the oscillating potential $mg_{E}z\,\lambda\sin\omega t$ and the
boson-boson interactions. They can be written in the form

$\displaystyle\frac{\partial}{\partial t}\widetilde{\eta}_{k}$
$\displaystyle=-i\epsilon_{k}\widetilde{\eta}_{k}-i\frac{g}{\hbar}{\displaystyle\sum\limits_{n}}E_{k,n}\,\widetilde{\eta}_{n}-i\frac{g}{\hbar}{\displaystyle\sum\limits_{n}}F_{k,n}\,\widetilde{\eta}_{n}$
(52) $\displaystyle\frac{\partial}{\partial t}\widetilde{\eta}_{k}^{+}$
$\displaystyle=+i\epsilon_{k}\widetilde{\eta}_{k}^{+}+i\frac{g}{\hbar}{\displaystyle\sum\limits_{n}}E_{k,n}^{+}\,\widetilde{\eta}_{n}^{+}+i\frac{g}{\hbar}{\displaystyle\sum\limits_{n}}F_{k,n}^{\ast}\,\widetilde{\eta}_{n}^{+}$
(53)

where

$\displaystyle E_{k,n}$
$\displaystyle={\displaystyle\int}dz\,\xi_{k}^{\ast}(z)\left(\widetilde{\psi}^{+}(z,t)\,\widetilde{\psi}(z,t)\,-\delta_{C}(z,z)\right)\,\xi_{n}(z)$
$\displaystyle E_{k,n}^{+}$
$\displaystyle={\displaystyle\int}dz\,\xi_{k}(z)\left(\widetilde{\psi}^{+}(z,t)\,\widetilde{\psi}(z,t)\,-\delta_{C}(z,z)\right)\,\xi_{n}^{\ast}(z)$
(54) $\displaystyle F_{k,n}$
$\displaystyle={\displaystyle\int}dz\,\xi_{k}^{\ast}(z)\left(mg_{e}z\,\lambda\sin\omega
t\right)\,\xi_{n}(z)$ (55)

Note that $F_{k,n}$ is periodic with period $T$. The time dependence of the
$E_{k,n}$ and $E_{k,n}^{+}$ will arise from the time dependence of the
$\widetilde{\eta}_{k}$ and $\widetilde{\eta}_{k}^{+}$ via the stochastic field
functions, and their time dependence will ultimately depend on $T$ and the
$\epsilon_{k}$. In this method the stochastic fields are determined at each
time point from Eq. (48), which then can also be used to determine the quantum
depletion (see Eq. (38)).

##### Initial Conditions - Gravitational Modes

Initial conditions in terms of gravitational modes are

$\widetilde{\psi}(z,0)={\displaystyle\sum\limits_{k}}\widetilde{\eta}_{k}(0)\,\xi_{k}(z)\qquad\widetilde{\psi}^{+}(z,0)={\displaystyle\sum\limits_{k}}\widetilde{\eta}_{k}^{+}(0)\,\xi_{k}^{\ast}(z)$
(56)

This requires a consideration of how to treat the preparation process, as will
be explained below.

####  2.3.3 Wannier Modes

Yet another choice for orthonormal mode functions that could be used to treat
the time crystal topic is the Wannier modes. These are defined as linear
combinations of the Floquet modes with frequencies in the first Floquet zone,
and are analogous to the Wannier functions defined for spatially periodic
potentials as linear combinations of Bloch functions associated with the first
Brillouin zone. Spatial Wannier functions are spatially localised around
different spatial lattice points; temporal Wannier functions are localised in
time around different time lattice points $nT$ (see Ref [30]). The Wannier
modes are

$\Phi_{nT}(z,t)=\mathcal{N}{\textstyle\sum\limits_{\\{\nu_{k}\\}}}\,\exp(-i\nu_{k}(t-nT))\;\phi_{k}(z,t)$
(57)

where the sum is over Floquet modes in the first Floquet zone and $N$ is a
normalizing factor, and $n$ is an integer. If $m$ is another integer the
periodic properties of the Floquet modes leads to the result
$\Phi_{(n+m)T}(z,t+mT)=\Phi_{nT}(z,t)$, which shows that for a given position
$z$ the Wannier mode $\Phi_{nT}(z,t)$ is a function of $(t-nT)$. This
indicates that the Wannier mode is centred in time around $nT$. It may be
temporally localised.

The Wannier modes satisfy the approximate orthogonality condition

${\textstyle\int}dz\,\Phi_{nT}(z,t)^{\ast}\,\Phi_{mT}(z,t)=\mathcal{N}^{2}{\textstyle\sum\limits_{\\{\nu_{k}\\}}}\exp(-i\nu_{k}(m-n)T)\approx\delta_{n,m}$
(58)

if the first Floquet zone is broad.

Approximate versions of Wannier modes are also used (Ref. [1]), based on a
two-mode approximation to the solutions of the Gross-Pitaevski equation. They
are two linear combinations of two Floquet modes designated as $\phi_{1}(z,t)$
and $\phi_{2}(z,t)$. For effects involving a period $sT$ we choose

$\displaystyle\Phi_{1}(z,t)$
$\displaystyle=\frac{1}{\sqrt{2}}(\phi_{1}(z,t)+\exp(-i\frac{2\pi}{sT}t)\,\phi_{2}(z,t))$
$\displaystyle\Phi_{2}(z,t)$
$\displaystyle=\frac{1}{\sqrt{2}}(\phi_{1}(z,t)-\exp(-i\frac{2\pi}{sT}t)\,\phi_{2}(z,t))$
(59)

These functions repeat over a period $sT$ and are orthogonal. The two Floquet
modes may be chosen to be similar to the condensate mode function at $t=0,$
enabling the initial conditions to be described in terms of these approximate
Wannier modes plus non-condensate modes such as Bogoliubov modes.

###  2.4 Preparation of the BEC - Condensate and Bogoliubov Modes

As described in Sect. 1, the BEC is prepared in the standard way in a trap
potential before it is allowed to bounce on the oscillating mirror. The
description of the states for the initially prepared BEC can be treated as a
time-independent problem via variational methods based on minimising the mean
value $\left\langle\widehat{H}\right\rangle$ of the energy, subject to
constraints such as the mean boson number
$\left\langle\widehat{N}\right\rangle$ being $N$. Here, we are only interested
in being able to represent the initial BEC state in terms of the Wigner
distribution just after the trap is switched off. The initial state of the BEC
for which $N_{c}$ bosonic atoms mainly occupy a single condensate mode may be
described theoretically in several ways, including the single-mode
approximation. In the single-mode approximation all the bosons are assumed to
occupy just one mode. Here the condensate bosons are described via a mean-
field approach, with the condensate wave function $\Phi_{c}(z)$ obtained from
the variational approach as the solution of a time-independent GPE, but with
the potential energy term now given by the time-independent trapping potential
$V_{trap}(z)$ and including the chemical potential $\mu$ to allow for the
constraint that $\left\langle\widehat{N}\right\rangle=N$ (see Eqs (60), (89)
below). However, collisional interactions cause bosons to be transferred from
the condensate mode to non-condensate modes, so the treatment of the initial
state should allow for this.

One widely used approach is to use Bogoliubov modes [31], [32], [25], [22],
[24], [26], [27] to describe the non-condensate modes. The field operator is
written as the sum of condensate and fluctuation fields. For the fluctuation
field Bogoliubov mode operators and their mode functions are introduced to
enable an approximation to the Hamiltonian in Eq (1) that is correct to
second-order in the fluctuation field to be expressed as the sum of a
condensate mode term and terms describing non-interacting quantum harmonic
oscillators (see Appendix 10 for details). First-order terms in the
fluctuation field are eliminated since the condensate field satisfies the
time-independent GPE. The mode functions are determined from Bogoliubov- de-
Gennes (BDG) equations based on the condensate wave function $\Phi_{c}$.
Bogoliubov modes thus allow for quantum fluctuations from the mean-field
theory described by the GPE and the condensate wave-function. This model will
be used to determine the initial stochastic quantities. The description in
terms of Bogoliubov modes for the BEC preparation is then matched to the
initial conditions for the BEC falling onto the oscillating mirror, where the
periodically driven BEC is treated as above in terms of Floquet or
gravitational modes.

####  2.4.1 Condensate Wave Function, Condensate and Non-Condensate Fields

The condensate wave-function $\Phi_{c}(z)$ is determined from the time-
independent GPE

$\left[-\frac{\hslash^{2}}{2m}\frac{\partial^{2}}{\partial
z^{2}}+V_{trap}(z)+g\\{\Phi_{c}(z)^{\ast}\Phi_{c}(z)\\}\right]\,\Phi_{c}(z)=\mu\,\Phi_{c}(z)$
(60)

where $\mu$ is the chemical potential. The condensate wave-function is
normalised in terms of the number $N$ of bosons in the condensate, as
${\displaystyle\int}dz\,\Phi_{c}(z)^{\ast}\Phi_{c}(z)=N$.

In Bogoliubov theory the field operator $\hat{\Psi}(z)$ is first written as
the sum of the condensate field $\widehat{\Psi}_{c}(z)$ and a fluctuation (or
non-condensate) field $\delta\widehat{\Psi}(z)$ as

$\hat{\Psi}(z)=\widehat{\Psi}_{c}(z)+\delta\widehat{\Psi}(z)$ (61)

where the condensate field and the fluctuation field are given by

$\displaystyle\widehat{\Psi}_{c}(z)$
$\displaystyle=\widehat{c}_{0}\psi_{c}(z)$ (62)
$\displaystyle\delta\widehat{\Psi}(z)$ $\displaystyle=\sum_{i\neq
C}\widehat{c}_{i}\psi_{i}(z)$ (63)

The condensate mode annihilation operator is $\widehat{c}_{0}$ and the
condensate mode function is

$\psi_{c}(z)=\Phi_{c}(z)/\sqrt{N}$ (64)

which is normalised to unity [32]. The equation for the condensate mode
function is

$\left[-\frac{\hslash^{2}}{2m}\frac{\partial^{2}}{\partial
z^{2}}+V_{trap}(z)+gN\\{\psi_{c}(z)^{\ast}\psi_{c}(z)\\}\right]\,\psi_{c}(z)=\mu\,\psi_{c}(z)$
(65)

This equation is used to eliminate terms in the Hamiltonian that are linear in
the fluctuation field. Note that the time variable $t$ is no longer present.

The fluctuation field may be expanded in terms of any suitable standard set of
normalised mode functions $\psi_{i}(z)$, where the corresponding non-
condensate mode annihilation operator is $\widehat{c}_{i}$. The non-condensate
modes are required to be orthogonal to $\psi_{c}(z)$, so we have the
constraint

$\int dz\,\psi_{c}(z)^{\ast}\delta\widehat{\Psi}(z)=0$ (66)

It follows that the non-zero commutation result for the fluctuation field is

$[\delta\widehat{\Psi}(z),\delta\widehat{\Psi}(z^{\\#})^{{\dagger}}]=\delta(z-z^{\\#})-\psi_{c}(z)\psi_{c}(z^{\\#})^{\ast}$
(67)

The choice of the non-condensate mode functions is arbitrary so far, but one
choice would be the normalised eigenfunctions of the _Gross-Pitaevskii
operator_ $\left(-\frac{\hslash^{2}}{2m}\frac{\partial^{2}}{\partial
z^{2}}+V_{trap}(z)+g\\{\Phi_{c}(z)^{\ast}\Phi_{c}(z)\\}\right)$ with
eigenvalues $\mu_{i}\neq\mu$. However, in the present paper we will use
Bogoliubov modes $u_{k}(z),v_{k}(z)$ rather than Gross-Pitaevski modes.

####  2.4.2 Bogoliubov Modes

The fluctuation field operator is expanded in terms of annihilation, creation
operators $\widehat{b}_{k},\widehat{b}_{k}^{{\dagger}}$ for Bogoliubov modes
and associated mode functions $u_{k}(z),v_{k}(z)$ as

$\delta\widehat{\Psi}(z)={\displaystyle\sum\limits_{k\neq
0}}\left[u_{k}(z)\,\widehat{b}_{k}-v_{k}(z)^{\ast}\,\widehat{b}_{k}^{{\dagger}}\right]\qquad\delta\widehat{\Psi}(z)^{{\dagger}}={\displaystyle\sum\limits_{k\neq
0}}\left[-v_{k}(z)\,\widehat{b}_{k}+u_{k}(z)^{\ast}\,\widehat{b}_{k}^{{\dagger}}\right]$
(68)

The summation over $k\neq 0$ is to indicate that the condensate mode
$\widehat{c}_{0}$ is excluded, since the fluctuation field is only intended to
treat non-condensate modes. The orthogonality condition (66) then leads to the
following orthogonality conditions between $u_{k}(z),v_{k}(z)$ and the
condensate mode function.

$\displaystyle\int dz\,\psi_{c}(z)^{\ast}u_{k}(z)$ $\displaystyle=0$
$\displaystyle\int dz\,\psi_{c}(z)\,v_{k}(z)$ $\displaystyle=0$ (69)

The commutation rules (67) for
$\delta\widehat{\Psi}(z),\delta\widehat{\Psi}(z^{\\#})^{{\dagger}}$ together
with the requirement that $\widehat{b}_{k},\widehat{b}_{k}^{{\dagger}}$
satisfy standard Bose commutation rules lead to the following conditions for
the $u_{k}(z),v_{k}(z)$

$\displaystyle{\displaystyle\sum\limits_{k\neq
0}}\left[u_{k}(z)\,v_{k}(z^{\\#})^{\ast}-v_{k}(z)^{\ast}\,u_{k}(z^{\\#})\right]$
$\displaystyle=0$ $\displaystyle{\displaystyle\sum\limits_{k\neq
0}}\left[u_{k}(z)\,u_{k}(z^{\\#})^{\ast}-v_{k}(z)^{\ast}\,v_{k}(z^{\\#})\right]$
$\displaystyle=\delta(z-z^{\\#})-\psi_{c}(z)\psi_{c}(z^{\\#})^{\ast}$ (70)

The inverse relations giving the $\widehat{b}_{k},\widehat{b}_{k}^{{\dagger}}$
in terms of the fluctuation fields are

$\displaystyle\widehat{b}_{k}$
$\displaystyle={\displaystyle\int}dz\,\left[u_{k}(z)\,\delta\widehat{\Psi}(z)+v_{k}(z)^{\ast}\,\delta\widehat{\Psi}(z)^{{\dagger}}\right]\qquad
k\neq 0$ $\displaystyle\widehat{b}_{k}^{{\dagger}}$
$\displaystyle={\displaystyle\int}dz\,\left[v_{k}(z)\,\delta\widehat{\Psi}(z)+u_{k}(z)^{\ast}\,\delta\widehat{\Psi}(z)^{{\dagger}}\right]\qquad
k\neq 0$ (71)

which is shown using Eqs (68), (70) and (66). Further conditions on
$u_{k}(z),v_{k}(z)$ follow using the commutation rules for
$\widehat{b}_{k},\widehat{b}_{k}^{{\dagger}}$ and
$\delta\widehat{\Psi}(z),\delta\widehat{\Psi}(z^{\\#})^{{\dagger}}$ and the
orthogonality conditions (69) in conjunction with Eq. (71). These are

$\displaystyle{\displaystyle\int}dz\,\left[u_{k}(z)\,v_{l}(z)\,-v_{k}(z)\,u_{l}(z)\right]$
$\displaystyle=0\qquad k,l\neq 0$
$\displaystyle{\displaystyle\int}dz\,\left[u_{k}(z)^{\ast}\,u_{l}(z)\,-v_{k}(z)^{\ast}\,v_{l}(z)\right]$
$\displaystyle=\delta_{k,l}\qquad k,l\neq 0$ (72)

This shows that the Bogoliubov mode functions $u_{k}(z),v_{k}(z)$ do not
satisfy standard orthogonality and normalisation conditions -
_biorthogonality_ conditions apply instead [26]. For each Bogoliubov frequency
$\omega_{k}$ we note that there are two mode functions involved.

####  2.4.3 Generalised Bogoliubov-de Gennes Equations

The Bogoliubov mode functions $u_{k}(z),v_{k}(z)$ and the mode frequencies
$\omega_{k}$ are chosen to satisfy a _generalised_ form [25] of the
Bogoliubov-de Gennes (BDG) equations

$\displaystyle(\mathcal{L}+g\,n_{C}(z))\,u_{k}(z)-g\,\Phi_{c}(z)^{2}\,v_{k}(z)$
$\displaystyle=\hslash\omega_{k}\,u_{k}(z)+C_{k}\psi_{c}(z)$
$\displaystyle-g\,\Phi_{c}^{\ast}(z)^{2}\,u_{k}(z)+(\mathcal{L}+g\,n_{C}(z))\,v_{k}(z)$
$\displaystyle=-\hslash\omega_{k}\,v_{k}(z)-C_{k}\psi_{c}^{\ast}(z)$ (73)

where $k\neq 0$. The differential operator $\mathcal{L}$ is defined by

$\mathcal{L}=\left[-\frac{\hslash^{2}}{2m}\frac{\partial^{2}}{\partial
z^{2}}+V_{trap}(z)+g\,n_{C}(z)-\mu\right]$ (74)

and the quantity $C_{k}$ is given by

$C_{k}=C_{k}[u_{k}(z),v_{k}(z)]=\int
dz\,g\,n_{C}(z)\,(\psi_{c}^{\ast}(z)\,u_{k}(z)-\psi_{c}(z)\,v_{k}(z))$ (75)

Note that the GPE (60) can be written
$\mathcal{L}\,\Phi_{c}(z)=\mathcal{L}\,\psi_{c}(z)=0$. The generalised BDG
equations are integro-differential equations, since $C_{k}$ is a functional of
$u_{k}(z),v_{k}(z)$. The equations depend only on quantities obtained from the
condensate wave function $\Phi_{c}(z)$, such as $n_{C}(z)$, $\Phi_{c}(z)^{2}$
and $\psi_{c}(z)$. Fortunately, the $u_{k}(z),v_{k}(z)$ and the Bogoliubov
frequencies $\omega_{k}$ can be obtained from eigenmodes of the standard BDG
equations, which are eigenvalue equations only involving differential
operators. Note that the sign convention differs from that in Ref [25].

It can be shown from the generalised BDG equations (73) that the eigenvalues
$\omega_{k}$ are real, and that the $u_{k}(z),v_{k}(z)$ satisfy the required
orthogonality conditions (69) with the condensate mode function $\psi_{c}(z)$.
The inclusion of the term involving $C_{k}$ is necessary to ensure
orthogonality. By convention $\omega_{k\text{ }}$is taken to be positive. The
proofs make use of the Hermitian properties of $\mathcal{L}$, where
${\displaystyle\int}dz\,U^{\ast}\,(\mathcal{L}\,V)={\displaystyle\int}dz\,(\mathcal{L}\,U^{\ast})\,V$
and $U,V$ are arbitrary functions. Furthermore, if $u_{k},v_{k}$ are solutions
for $\omega_{k}$, then $v_{k}^{\ast},u_{k}^{\ast}$ are solutions for
$-\omega_{k}$. The proof involves the property of $C_{k}$ that
$C_{k}[u_{k}(z),v_{k}(z)]^{\ast}=-C_{k}[v_{k}^{\ast}(z),u_{k}^{\ast}(z)]$. The
generalised BDG equations seem to have a solution for $\omega_{0}=0$ with
$u_{0}(z)=\psi_{c}(z)$ and $v_{0}(z)=\psi_{c}^{\ast}(z)$. However, this is
inconsistent with Eq. (69), based on the non-condensate field being required
to be orthogonal to the condensate mode (see Eq.(66)).

The second-order approximation to the Hamiltonian (1) can be written in terms
of the Bogoliubov mode operators $\widehat{b}_{k},\widehat{b}_{k}^{{\dagger}}$
as the sum of uncoupled harmonic oscillators, one for each Bogoliubov mode.
The form is given in Appendix 10 as Eq.(136).

####  2.4.4 Standard Bogoliubov-de Gennes Equations

The standard Bogoliubov-de Gennes equations can be obtained by replacing the
modes $u_{k},v_{k}$ by new modes $\overline{u}_{k},\overline{v}_{k}$ via the
expressions

$u_{k}=\overline{u}_{k}-\frac{C_{k}\psi_{c}}{\hslash\omega_{k}}\qquad\qquad
v_{k}=\overline{v}_{k}-\frac{C_{k}\psi_{c}^{\ast}}{\hslash\omega_{k}}$ (76)

Substituting into the generalised BDG equations we then obtain the _standard_
BDG equations [21]

$\displaystyle(\mathcal{L}+g\,n_{C}(z))\,\overline{u}_{k}(z)-g\,\Phi_{c}(z)^{2}\,\overline{v}_{k}(z)$
$\displaystyle=\hslash\omega_{k}\,\overline{u}_{k}(z)$
$\displaystyle-g\,\Phi_{c}^{\ast}(z)^{2}\,\overline{u}_{k}(z)+(\mathcal{L}+g\,n_{C}(z))\,\overline{v}_{k}(z)$
$\displaystyle=-\hslash\omega_{k}\,\overline{v}_{k}(z)$ (77)

where the orthogonality conditions (69) have been used. The quantity $C_{k}$
is now also given by

$C_{k}=\int
dz\,g\,n_{C}(z)\,(\psi_{c}^{\ast}(z)\,\overline{u}_{k}(z)-\psi_{c}(z)\,\overline{v}_{k}(z))$
(78)

in terms of these modes. A generalised orthogonality condition $\int
dz\,(\psi_{c}^{\ast}(z)\,\overline{u}_{k}(z)-\psi_{c}(z)\,\overline{v}_{k}(z))=0$
follows from the standard BDG equations, but this does not mean that the new
modes $\overline{u}_{k},\overline{v}_{k}$ individually satisfy the required
orthogonality conditions (69) with the condensate mode function.

We can therefore obtain the true $u_{k},v_{k}$ by using the standard BDG
equations (77) to first determine the $\overline{u}_{k},\overline{v}_{k}$ and
then use (76) to determine the $u_{k},v_{k}$. Since the $u_{k},v_{k}$ satisfy
the orthogonality conditions (69), we can use this feature to determine the
quantity $C_{k}/\hslash\omega_{k}$. We find that

$D_{k}=\frac{C_{k}}{\hslash\omega_{k}}=\int
dz\,\psi_{c}^{\ast}(z)\,\overline{u}_{k}(z)=\int
dz\,\psi_{c}(z)\,\overline{v}_{k}(z)$ (79)

####  2.4.5 Field Operators - Condensate and Bogoliubov Modes

This gives the following expressions for the field operators after introducing
the $\overline{u}_{k},\overline{v}_{k}$ into the expressions (61) for the
fluctuation field. We have

$\displaystyle\hat{\Psi}(z)$
$\displaystyle=\widehat{c}_{0}\psi_{c}(z)+{\displaystyle\sum\limits_{k\neq
0}}\left[\left(\overline{u}_{k}(z)-D_{k}\psi_{c}(z)\right)\,\widehat{b}_{k}-\left(\overline{v}_{k}(z)-D_{k}\psi_{c}^{\ast}(z)\right)^{\ast}\,\widehat{b}_{k}^{{\dagger}}\right]$
$\displaystyle\hat{\Psi}(z)^{{\dagger}}$
$\displaystyle=\widehat{c}_{0}^{{\dagger}}\psi_{c}^{\ast}(z)+{\displaystyle\sum\limits_{k\neq
0}}\left[-\left(\overline{v}_{k}(z)-D_{k}\psi_{c}^{\ast}(z)\right)\,\widehat{b}_{k}+\left(\overline{u}_{k}(z)-D_{k}\psi_{c}(z)\right)^{\ast}\,\widehat{b}_{k}^{{\dagger}}\right]$
(80)

Similar equations to (80) are set out in Ref. [21]. We can use these
expressions to relate the stochastic field functions at the initial time. The
Bogoliubov modes $\overline{u}_{k},\overline{v}_{k}$ involved are those
obtained from the standard BDG equations (77).

####  2.4.6 Field Functions - Condensate and Bogoliubov Modes

As an alternative to Floquet or gravitational modes, the field functions
$\Psi(z)$, $\Psi^{+}(z)$ and the stochastic field functions can also be
expanded in terms of Bogoliubov modes. The mode annihilation, creation
operators $\widehat{b}_{k},\widehat{b}_{k}^{{\dagger}}$ are represented by
phase-space variables $\beta_{k},\beta_{k}^{+}$, which are replaced by
stochastic variables $\widetilde{\beta}_{k},\widetilde{\beta}_{k}^{+}$.
Similarly, the condensate mode annihilation, creation operators
$\widehat{c}_{0},\widehat{c}_{0}^{{\dagger}}$ are represented by phase-space
variables $\gamma_{0},\gamma_{0}^{+}$, which are replaced by stochastic
variables $\widetilde{\gamma}_{0},\widetilde{\gamma}_{0}^{+}$. Based on Eqs.
(80) and re-introducing the original Bogoliubov modes we have

$\displaystyle\psi(z)$
$\displaystyle=\gamma_{0}\psi_{c}(z)+{\displaystyle\sum\limits_{k\neq
0}}\left[u_{k}(z)\,\beta_{k}-v_{k}(z)^{\ast}\,\beta_{k}^{+}\right]$
$\displaystyle\psi^{+}(z)$
$\displaystyle=\gamma_{0}^{+}\psi_{c}^{\ast}(z)+{\displaystyle\sum\limits_{k\neq
0}}\left[-v_{k}(z)\,\beta_{k}+u_{k}(z)^{\ast}\,\beta_{k}^{+}\right]$ (81)
$\displaystyle\widetilde{\psi}(z)$
$\displaystyle=\widetilde{\gamma}_{0}\psi_{c}(z)+{\displaystyle\sum\limits_{k\neq
0}}\left[u_{k}(z)\,\widetilde{\beta}_{k}-v_{k}(z)^{\ast}\,\widetilde{\beta}_{k}^{+}\right]$
$\displaystyle\widetilde{\psi}^{+}(z)$
$\displaystyle=\widetilde{\gamma}_{0}^{+}\psi_{c}^{\ast}(z)+{\displaystyle\sum\limits_{k\neq
0}}\left[-v_{k}(z)\,\widetilde{\beta}_{k}+u_{k}(z)^{\ast}\,\widetilde{\beta}_{k}^{+}\right]$
(82)

where $u_{k}(z),v_{k}(z)$ are given in Eqs, (76) and (79). Similar equations
to (82) are set out in Ref. [21]. Note here that none of the stochastic
quantities are time dependent.

###  2.5 Stochastic Field Functions during Driven Evolution

During the evolution of the system of interacting bosons driven by the
oscillating potential and bouncing under the influence of the gravitation
field, the stochastic field functions and their expansion coefficients (or
stochastic mode amplitudes) are time dependent. In the cases of Floquet modes
and gravitational modes the stochastic field functions are given in Eqs. (45)
and (48), respectively.

However, during evolution the field operators, field functions and stochastic
field functions can also be expressed in terms of the time-dependent
condensate mode function $\psi_{c}(z,t)$ and any set of orthogonal non-
condensate mode functions $\psi_{i}(z,t)$ \- as in Eq. (63) for the
fluctuation field operator. Thus, for the stochastic field functions we have
(see Eq. (34))

$\displaystyle\widetilde{\psi}(z,t)$
$\displaystyle=\widetilde{\gamma}_{0}(t)\psi_{c}(z,t)+{\displaystyle\sum\limits_{i\neq
0}}\widetilde{\gamma}_{i}(t)\psi_{i}(z,t)$
$\displaystyle\widetilde{\psi}^{+}(z,t)$
$\displaystyle=\widetilde{\gamma}_{0}^{+}(t)\psi_{c}^{\ast}(z,t)+{\displaystyle\sum\limits_{i\neq
0}}\widetilde{\gamma}_{i}^{+}(t)\psi_{i}^{\ast}(z,t)$ (83)

The stochastic mode amplitudes are
$\widetilde{\gamma}_{0}(t),\widetilde{\gamma}_{0}^{+}(t)$ for the condensate
mode and $\widetilde{\gamma}_{i}(t),\widetilde{\gamma}_{i}^{+}(t)$ for the
non-condensate modes. These expressions can be equated to those in Eqs. (82),
(45) and (48) at $t=0$. In regard to the last form we have
$\widetilde{\gamma}_{0}(0)=\widetilde{\gamma}_{0},$
$\widetilde{\gamma}_{0}^{+}(0)=\widetilde{\gamma}_{0}^{+}$ in view of
$\psi_{c}(z,0)=\psi_{c}(z)$ \- hence the same notation.

####  2.5.1 Initial Conditions for Stochastic Amplitudes of Gravitational
Modes

The initial conditions for the gravitational mode approach can be obtained
from the Bogoliubov mode approach via equating Eqs. (48) and (82) at $t=0$ and
using the orthogonality properties of the gravitational modes. We find that

$\displaystyle\widetilde{\eta}_{k}(0)={\textstyle\int}dz\,\xi_{k}^{\ast}(z)\,\left[\widetilde{\gamma}_{0}\psi_{c}(z)+{\displaystyle\sum\limits_{k\neq
0}}\left[u_{k}(z)\,\widetilde{\beta}_{k}-v_{k}(z)^{\ast}\,\widetilde{\beta}_{k}^{+}\right]\right]$
$\displaystyle\widetilde{\eta}_{k}^{+}(0)={\textstyle\int}dz\,\xi_{k}(z)\,\left[\widetilde{\gamma}_{0}^{+}\psi_{c}^{\ast}(z)+{\displaystyle\sum\limits_{k\neq
0}}\left[-v_{k}(z)\,\widetilde{\beta}_{k}+u_{k}(z)^{\ast}\,\widetilde{\beta}_{k}^{+}\right]\right]$

where $u_{k}(z),v_{k}(z)$ are given in Eqs. (76) and (79). This shows how a
particular choice of stochastic phase-space variables
$\widetilde{\gamma}_{0},\widetilde{\gamma}_{0}^{+},$
$\widetilde{\beta}_{q},\widetilde{\beta}_{q}^{+}$ for the Bogoliubov mode
treatment of the BEC preparation can be translated into the choice of
stochastic phase-space variables
$\widetilde{\eta}_{k}(0),\widetilde{\eta}_{k}^{+}(0)$ for the gravitational
mode treatment of the periodically driven BEC. Essentially, the
$\widetilde{\eta}_{k}(0),\widetilde{\eta}_{k}^{+}(0)$ for the gravitational
mode treatment are linearly dependent on the
$\widetilde{\gamma}_{0},\widetilde{\gamma}_{0}^{+},\widetilde{\beta}_{q},\widetilde{\beta}_{q}^{+}$.

We then see that we can write the initial stochastic fields for gravitational
mode evolution as

$\widetilde{\psi}(z,0)=\widetilde{\gamma}_{0}(0)\psi_{c}(z)+\delta\widetilde{\psi}(z,0)\qquad\widetilde{\psi}^{+}(z,0)=\widetilde{\gamma}_{0}^{+}(0)\psi_{c}^{\ast}(z)+\delta\widetilde{\psi}^{+}(z,0)$
(85)

with

$\displaystyle\delta\widetilde{\psi}(z,0)$
$\displaystyle={\textstyle\sum\limits_{k}}\xi_{k}(z)\,{\textstyle\int}dz^{\\#}\,\xi_{k}^{\ast}(z^{\\#})\,\times{\displaystyle\sum\limits_{k\neq
0}}\left[u_{k}(z^{\\#})\,\widetilde{\beta}_{k}-v_{k}(z^{\\#})^{\ast}\,\widetilde{\beta}_{k}^{+}\right]$
$\displaystyle\delta\widetilde{\psi}^{+}(z,0)$
$\displaystyle={\textstyle\sum\limits_{k}}\xi_{k}^{\ast}(z)\,{\textstyle\int}dz^{\\#}\,\xi_{k}(z^{\\#})\,\times{\displaystyle\sum\limits_{k\neq
0}}\left[-v_{k}(z^{\\#})\,\widetilde{\beta}_{k}+u_{k}(z^{\\#})^{\ast}\,\widetilde{\beta}_{k}^{+}\right]$
(86)

where $u_{k}(z),v_{k}(z)$ are given in Eqs. (76) and (79). These expressions
allow for quantum fluctuations in both the condensate and non-condensate
modes.

Analogous expressions to (LABEL:Eq._StochFloqBogolRelns) apply for the initial
stochastic amplitudes $\widetilde{\alpha}_{k}(0)$,
$\widetilde{\alpha}_{k}^{+}(0)$ of Floquet modes with the replacement of the
gravitational modes $\,\xi_{k}(z)\,$by the Floquet modes $\phi_{k}(z,0)$, see
Appendix 9.

###  2.6 Quantum State for Initial BEC

The question now arises - Having described the prepared BEC via the condensate
mode and the Bogoliubov modes, what is the quantum density operator? This
density operator is needed to define the Wigner distribution function for this
state and hence the stochastic properties of the
$\widetilde{\beta}_{q}(0),\widetilde{\beta}_{q}^{+}(0)$ and
$\widetilde{\gamma}_{0}(0),\widetilde{\gamma}_{0}^{+}(0)$, which are
equivalent to this initial Wigner function. If the Bogoliubov modes are in
thermal equilibrium at temperature $T_{0}$ then the density operator for the
non-condensate modes should be given by

$\widehat{\rho}_{NC}=\exp(-{\textstyle\sum\limits_{k}}\hslash\gamma_{k}\,\widehat{b}_{k}^{{\dagger}}\widehat{b}_{k}/k_{B}T_{0})/Tr\left[\exp(-{\textstyle\sum\limits_{k}}\hslash\gamma_{k}\,\widehat{b}_{k}^{{\dagger}}\widehat{b}_{k}/k_{B}T_{0})\right]$.
This of course ignores interactions between the condensate mode and the non-
condensate Bogoliubov modes.

At $T_{0}=0$ the non-condensate Bogoliubov modes are unoccupied and this would
be approximately the case at temperatures where the BEC exists. Consistent
with the GPE that has been used to describe the condensate bosons, all $N_{c}$
bosons are assumed to be in the condensate mode - whose mode function is
obtained from the time-independent Gross-Pitaevskii equation (60) - and the
state involved is a Fock state. The combined state for both condensate and
non-condensate modes will then be a pure state whose density operator commutes
with the number operator, and hence complies with the requirement of being
phase invariant. Descriptions in which the condensate bosons are in a Glauber
coherent state are often used for mathematical convenience, but as such states
are not phase invariant they are unphysical. For completeness, this approach
is set out in Appendix 11.

In this section the time value is $t=0$, which is left understood for
simplicity.

####  2.6.1 Density Operator for Initial State

The density operator is given by

$\widehat{\rho}(0)=\left|\Phi_{N0}\right\rangle\left\langle\Phi_{N0}\right|$
(87)

where

$\displaystyle\left|\Phi_{N0}\right\rangle$
$\displaystyle=\left|N_{c}\right\rangle_{0}\times\left|0\right\rangle_{B}$
$\displaystyle\left|N\right\rangle_{0}$
$\displaystyle=\frac{(\widehat{c}_{0}^{{\dagger}})^{N}}{\sqrt{N!}}\left|0\right\rangle_{c}\qquad\left|0\right\rangle_{B}={\textstyle\prod\limits_{k\neq
0}}\,\left|0\right\rangle_{k}$ (88)

This state is the product of vacuum states for the Bogoliubov modes and a Fock
state with $N$ bosons for the condensate mode.

####  2.6.2 First and Second-Order QCF for Initial State

We then find that

$\displaystyle\hat{\Psi}(z)\,\left|\Phi_{N0}\right\rangle$
$\displaystyle=\sqrt{N}\,\psi_{c}(z)\,\left|N-1\right\rangle_{0}\times\left|0\right\rangle_{B}+\left|N\right\rangle_{0}\times{\displaystyle\sum\limits_{k\neq
0}}(-v_{k}^{\ast}(z)\,\left|0\right\rangle_{1}..\left|1\right\rangle_{k}...\left|0\right\rangle_{n}$
$\displaystyle\left\langle\hat{\Psi}(z)\right\rangle$ $\displaystyle=0$
$\displaystyle\left\langle\hat{\Psi}(z)^{{\dagger}}\,\hat{\Psi}(z)\right\rangle$
$\displaystyle=n_{C}(z)=N\,\psi_{c}(z)^{\ast}\psi_{c}(z)+{\displaystyle\sum\limits_{k\neq
0}}(v_{k}^{\ast}(z)\,v_{k}(z))\,\approx\Phi_{c}(z)^{\ast}\Phi_{c}(z)$
$\displaystyle\left\langle\hat{\Psi}(z^{\\#})^{{\dagger}}\,\hat{\Psi}(z)\right\rangle$
$\displaystyle=P(z,z^{\\#},0)=N\,\psi_{c}(z^{\\#})^{\ast}\psi_{c}(z)+{\displaystyle\sum\limits_{k\neq
0}}(v_{k}^{\ast}(z)\,v_{k}(z^{\\#}))\,\approx\Phi_{c}(z^{\\#})^{\ast}\Phi_{c}(z)$
(89)

where the terms involving $v_{k}(z)$ can be neglected since $N\gg 1$. We see
that the mean value of the field operator is zero, as required for a phase-
invariant state. However, the mean value of the number density operator is
still obtained from the condensate wave-function. The Fock state is a good
description of the initial BEC state, and will therefore be used to determine
the initial conditions. Note that the last equation demonstrates long-range
spatial order in the BEC, since the QCF will be non-zero for $|z-z^{\\#}|$
within the extent of the condensate wave function. The result
$P(z,z^{\\#},0)=N\,\psi_{c}(z^{\\#})^{\ast}\psi_{c}(z)$ expresses the QCF in
terms of natural orbitals - in this case the only natural orbital occupied is
the condensate mode function, where the occupancy is $N_{c}$.

####  2.6.3 Wigner Distribution Functions for Initial State

For any given initial state, the stochastic averages of products of stochastic
field functions $\widetilde{\psi}(z),\widetilde{\psi}^{+}(z)$ based on phase-
space theory using the Wigner distribution functional can be obtained from the
quantum correlation functions for normally ordered products of field operators
$\hat{\Psi}(z),\hat{\Psi}(z)^{{\dagger}}$ and then using Wick’s theorem [33]
to obtain the quantum correlation functions for symmetrically ordered products
of the same field operators. These QCF for the symmetrically ordered products
are then equal to the stochastic averages of the products of stochastic field
functions associated with the original field operators. From these QCF the
stochastic averages of products of stochastic phase-space variables
$\widetilde{\gamma}_{0}(0),\widetilde{\gamma}_{0}^{+}(0),$
$\widetilde{\beta}_{q}(0),\widetilde{\beta}_{q}^{+}(0)$ for the Bogoliubov
mode treatment of the BEC preparation could be obtained. This method would be
required when the condensate mode and the non-condensate modes are correlated
and the initial density operator does not factorise into separate density
operators for the condensate and non-condensate modes. For completeness, the
more general method is described in Appendix 12.

However, for the initial state given by Eq. (87) the density operator
factorises into a density operator $\widehat{\rho}_{C}(0)$ for the condensate
mode and a density operator $\widehat{\rho}_{NC}(0)$ for the non-condensate
Bogoliubov modes, which are all in their vacuum states

$\displaystyle\widehat{\rho}(0)$
$\displaystyle=\widehat{\rho}_{C}(0)\otimes\widehat{\rho}_{NC}(0)$
$\displaystyle\widehat{\rho}_{C}(0)$
$\displaystyle=\left|\Phi_{N0}\right\rangle\left\langle\Phi_{N0}\right|\qquad\widehat{\rho}_{NC}(0)={\textstyle\prod\limits_{k\neq
0}}(\left|0\right\rangle\left\langle 0\right|)_{k}$ (90)

For this situation we can show that the Wigner distribution function
$W(\gamma_{0},\gamma_{0}^{+},\mathbf{\beta},\mathbf{\beta}^{+})$ factorises
into separate Wigner distribution functions $W_{C}(\gamma_{0},\gamma_{0}^{+})$
and $W_{NC}(\mathbf{\beta},\mathbf{\beta}^{+})$ for the condensate mode and
the non-condensate Bogoliubov modes, where
$\mathbf{\beta}\equiv\\{\beta_{1},...\beta_{k},..)$ and
$\mathbf{\beta}^{+}\equiv\\{\beta_{1}^{+},...\beta_{k}^{+},..)$ We can then
calculate the stochastic averages of products of stochastic phase-space
variables by considering the condensate and non-condensate modes separately.

The proof of the result

$W(\gamma_{0},\gamma_{0}^{+},\mathbf{\beta},\mathbf{\beta}^{+})=W_{C}(\gamma_{0},\gamma_{0}^{+})\times
W_{NC}(\mathbf{\beta},\mathbf{\beta}^{+})$ (91)

where the Wigner functions are normalised to unity -
${\textstyle\int}d^{2}\gamma_{0}d^{2}\gamma_{0}^{+}\,W_{C}(\gamma_{0},\gamma_{0}^{+})=1$
and
${\textstyle\int}d^{2}\mathbf{\beta}\,d^{2}\mathbf{\beta}^{+}\,W_{NC}(\mathbf{\beta},\mathbf{\beta}^{+})=1$
\- is given in Appendix 13. This result means that phase-space averages that
give QCF for products of condensate and non-condensate mode operators will
factorise into a separate phase-space average for the condensate mode and a
phase-space average for the non-condensate modes. This in turn means that the
equivalent stochastic average of products of stochastic phase-space variables
for the condensate mode and stochastic phase-space variables for the non-
condensate modes will factorise into separate stochastic averages for the
condensate and non-condensate stochastic variables. Thus, for example,
$\overline{\widetilde{\gamma}_{0}\widetilde{\beta}_{k}^{+}}=\overline{\widetilde{\gamma}_{0}}$
$\times\overline{\widetilde{\beta}_{k}}$ .

The result uses Morgan’s approach [25] and involves first establishing a
relation between the standard non-condensate mode operators
$\widehat{c}_{i},\widehat{c}_{i}^{{\dagger}}$ that were introduced in Eq. (63)
and the Bogoliubov mode operators
$\widehat{b}_{k},\widehat{b}_{k}^{{\dagger}}$ which are given in matrix form
(the convention differs from that in Ref [25]) as

$\left[\begin{array}[c]{c}\widehat{\mathbf{c}}\\\
\widehat{\mathbf{c}}^{{\dagger}}\end{array}\right]=\left[\begin{array}[c]{cc}U&-V\\\
-V^{\ast}&U^{\ast}\end{array}\right]\times\left[\begin{array}[c]{c}\widehat{\mathbf{b}}\\\
\widehat{\mathbf{b}}^{{\dagger}}\end{array}\right]$ (92)

where the matrices $U,V$ have elements given by

$U_{i,k}={\textstyle\int}dz\,\psi_{i}(z)^{\ast}\,u_{k}(z)\qquad
V_{i,k}={\textstyle\int}dz\,\psi_{i}(z)^{\ast}\,v_{k}(z)^{\ast}$ (93)

This is a linear canonical transformation in which commutation rules are
preserved. For the matrices, $T$ denotes transverse, $\ast$ denotes complex
conjugation and ${\dagger}$ denotes the Hermitian adjoint. The bosonic
commutation rules for both sets of non-condensate mode operators leads to the
following matrix equations.

$\displaystyle-U\,V^{T}+V\,U^{T}$ $\displaystyle=0$ $\displaystyle
U\,U^{{\dagger}}-V\,V^{{\dagger}}$ $\displaystyle=E$ (94)

Further features of these matrices are set out in Appendix 13.

####  2.6.4 Stochastic Averages - Condensate Mode

For the condensate mode the density operator is given by (87) and as
$\widehat{c}_{0}\left|N\right\rangle=\sqrt{N}\left|N-1\right\rangle$ it is
easy to show that the first and second-order QCF for normally ordered
condensate mode operators are given by

$\displaystyle\left\langle\widehat{c}_{0}\right\rangle$
$\displaystyle=\left\langle\widehat{c}_{0}^{{\dagger}}\right\rangle=0$ (95)
$\displaystyle\left\langle\widehat{c}_{0}\widehat{c}_{0}\right\rangle$
$\displaystyle=\left\langle\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}^{{\dagger}}\right\rangle=0\qquad\left\langle\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}\right\rangle=N$
(96)

Since
$\left\\{\widehat{c}_{0}\widehat{c}_{0}\right\\}=\widehat{c}_{0}\widehat{c}_{0}$,
$\left\\{\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}^{{\dagger}}\right\\}=\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}^{{\dagger}}$
and
$\left\\{\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}\right\\}=\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}+1/2$
we see that the stochastic averages are

$\displaystyle\overline{\widetilde{\gamma}_{0}}$
$\displaystyle=\overline{\widetilde{\gamma}_{0}^{+}}=0$ (97)
$\displaystyle\overline{\widetilde{\gamma}_{0}\widetilde{\gamma}_{0}}$
$\displaystyle=\overline{\widetilde{\gamma}_{0}^{+}\widetilde{\gamma}_{0}^{+}}=0\qquad\overline{\widetilde{\gamma}_{0}^{+}\widetilde{\gamma}_{0}}=N+\frac{1}{2}$
(98)

The distribution function generating the stochastic $\widetilde{\gamma}_{0}$
and $\widetilde{\gamma}_{0}^{+}$ thus must produce a mean value of zero for
the first order, and $N+\frac{1}{2}$ for the non-zero second-order case
$\overline{\widetilde{\gamma}_{0}^{+}\widetilde{\gamma}_{0}}$.

Higher order QCFs can also be considered. The third-order QCF for normally
ordered condensate mode operators are all zero. The only non-zero fourth-order
QCF for normally ordered condensate mode operators is
$\left\langle\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}\widehat{c}_{0}\right\rangle=N(N-1)$.
Since
$\\{\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}\widehat{c}_{0}\\}=(\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}\widehat{c}_{0}+\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}+\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}\widehat{c}_{0}\widehat{c}_{0}^{{\dagger}}+\widehat{c}_{0}\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}+\widehat{c}_{0}\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}\widehat{c}_{0}^{{\dagger}}+\widehat{c}_{0}\widehat{c}_{0}\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}^{{\dagger}})/6$
we can use Wick’s theorem [33] to show that
$\\{\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}\widehat{c}_{0}\\}=\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}\widehat{c}_{0}+2\widehat{c}_{0}^{{\dagger}}\widehat{c}_{0}+1/2$.
We then see that the non-zero stochastic average is

$\displaystyle\overline{(\widetilde{\gamma}_{0}^{+})^{2}(\widetilde{\gamma}_{0})^{2}}$
$\displaystyle=N(N-1)+2N+1/2$ $\displaystyle=N^{2}+N+1/2$ (99)

This result is given in Ref [34]. If we choose
$\widetilde{\gamma}_{0}=\sqrt{\widetilde{\Gamma}_{0}}\,\exp(i\widetilde{\phi}_{0})$
and
$\widetilde{\gamma}_{0}^{+}=\sqrt{\widetilde{\Gamma}_{0}}\,\exp(-i\widetilde{\phi}_{0})$,
where $\widetilde{\Gamma}_{0}$ is a normally distributed Gaussian random
variable with a mean equal to $N+\frac{1}{2}$ and a variance of $1/4$, and we
choose $\widetilde{\phi}_{0}$ to be a uniformly distributed random phase over
$0$ to $2\pi$, then all the results in Eqs.(97), (98) and (99) are obtained.
The Wigner distribution function for a Fock state is given in [34].

####  2.6.5 Stochastic Averages - Non-Condensate Modes

For the non-condensate mode the density operator is given by (90) and as
$\widehat{b}_{k}\left|0\right\rangle_{k}=0\times\left|0\right\rangle_{k}$ it
is easy to show that the first and second-order QCF for normally ordered non-
condensate mode operators are given by

$\displaystyle\left\langle\widehat{b}_{k}\right\rangle$
$\displaystyle=\left\langle\widehat{b}_{k}^{{\dagger}}\right\rangle=0$ (100)
$\displaystyle\left\langle\widehat{b}_{k}\widehat{b}_{l}\right\rangle$
$\displaystyle=\left\langle\widehat{b}_{k}^{{\dagger}}\widehat{b}_{l}^{{\dagger}}\right\rangle=0\qquad\left\langle\widehat{b}_{k}^{{\dagger}}\widehat{b}_{l}\right\rangle=0$
(101)

Since
$\left\\{\widehat{b}_{k}\widehat{b}_{l}\right\\}=\widehat{b}_{k}\widehat{b}_{l}$,
$\left\\{\widehat{b}_{k}^{{\dagger}}\widehat{b}_{l}^{{\dagger}}\right\\}=\widehat{b}_{k}^{{\dagger}}\widehat{b}_{l}^{{\dagger}}$
and
$\left\\{\widehat{b}_{k}^{{\dagger}}\widehat{b}_{l}\right\\}=\widehat{b}_{k}^{{\dagger}}\widehat{b}_{l}+(1/2)\,\delta_{k,l}$
we see that the stochastic averages are

$\displaystyle\overline{\widetilde{\beta}_{k}}$
$\displaystyle=\overline{\widetilde{\beta}_{k}^{+}}=0$ (102)
$\displaystyle\overline{\widetilde{\beta}_{k}\widetilde{\beta}_{l}}$
$\displaystyle=\overline{\widetilde{\beta}_{k}^{+}\widetilde{\beta}_{l}^{+}}=0\qquad\overline{\widetilde{\beta}_{k}^{+}\widetilde{\beta}_{l}}=\frac{1}{2}\,\delta_{k,l}$
(103)

The distribution function generating the stochastic $\widetilde{\beta}_{k}$
and $\widetilde{\beta}_{k}^{+}$ thus must produce a mean value of zero for the
first order, and $\frac{1}{2}$ for the non-zero second order case
$\overline{\widetilde{\beta}_{k}^{+}\widetilde{\beta}_{k}}$.

It is of some interest to calculate the stochastic averages for the stochastic
phase-space variables associated with the standard modes
$\widehat{c}_{i},\widehat{c}_{i}^{{\dagger}}$ associated with Eq. (63). It can
be shown (see Appendix 13) that the vacuum state for Bogoliubov modes is
equivalent to a squeezed vacuum state for the standard non-condensate modes.

##  3 Numerical Results

The numerics in this paper are based on the truncated Wigner approximation
(TWA) and to a lesser extent on the mean-field Gross-Pitaevskii equation
(GPE). Our aim in the present paper is to study the DTC with periodicity $sT$
for the case of period doubling $(s=2)$, where $T$ is the periodicity of the
oscillating mirror.

The position probability density (PPD)
$F(z,t)=Tr(\hat{\Psi}^{{\dagger}}(z)\hat{\Psi}(z)\rho(t))$ versus $z$ plots
for various times (Eq.(18)) depend either on the time-dependent stochastic
field functions $\widetilde{\psi}(z,t)$ in the TWA case or on the time-
dependent condensate wave function $\Phi_{c}(z,t)$, both considered as a
function of time (see Eqs.(8), (9), (10)). The same is true for the quantum
correlation function (QCF)
$P(z,z^{\\#},t)=Tr(\hat{\Psi}^{{\dagger}}(z^{\\#})\hat{\Psi}(z)\rho(t))$ (see
Eq.(20)) and the one-body projector (OBP) (both for the condensate mode
$M_{c}(t)={\displaystyle\int}{\displaystyle\int}dz\,dz^{\\#}\,\psi_{c}(z^{\\#},0)\,\,\psi_{c}^{\ast}(z,0)\,P(z,z^{\\#},t)$
and for the Wannier modes $\Phi_{1}$ and $\Phi_{2}$ (see Eq.(59), with $s=2$)
and its Fourier transform (FT) (see Eqs.(26), (28), (32)). The natural
orbitals and their occupation numbers are determined from the QCF (see
Eqs.(21), (22)). The mean energy is determined from the stochastic field
functions (see Eq. (13)) with the most important contributions given in
Eq.(51) using gravitational modes. The quantum depletion from the condensate
mode $N_{D}(t)=N-N_{C}(t)$ is determined from Eq. (41).

In both the TWA and GPE numerics, the behaviour of these functions is
deterministic and depends on the initial conditions $\widetilde{\psi}(z,0)$
and $\Phi_{c}(z,0)$. The difference between the TWA and GPE is that in the
former a stochastic ensemble of $\widetilde{\psi}(z,0)$ is considered whereas
in the latter just a single $\Phi_{c}(z,t)$ is involved (see Sect.2.2). In the
TWA numerics we find it convenient to choose stochastic field functions such
that $\widetilde{\psi}^{+}(z,t)=\widetilde{\psi}(z,t)^{\ast}$ and stochastic
phase amplitudes such that
$\widetilde{\eta}_{k}^{+}(t)=\widetilde{\eta}_{k}(t)^{\ast}$, etc.

For the initial conditions the BEC is assumed to be at zero temperature, with
all bosons in a single mode. For both the TWA and GPE cases the initial
condensate mode function is based on the time-independent GPE for a harmonic
trap with non-condensate modes being Bogoliubov modes. For the Wannier initial
condition a linear combination of two Wannier modes is taken as the initial
condensate mode function, with non-condensate modes being any orthogonal set
of modes. The linear combination is chosen to minimise the energy functional
as described in Ref. [1]. For the TWA case the initial condition allows for
unoccupied Bogoliubov modes in the harmonic trap case, and unoccupied non-
condensate modes in the Wannier case. There are no unoccupied non-condensate
modes in the GPE approach by definition. The loss of bosons from the
condensate mode is specified by the quantum depletion (see Eqs.(36), (37),
(38)). For the TWA case there is a stochastic ensemble of phase-space
amplitudes for the condensate and Bogoliubov modes (see Sect. 2.4.6 and Eq.
(82)). Equivalent initial stochastic phase-space amplitudes for gravitational
modes are determined by equating expressions for the initial stochastic field
functions (see Sect.2.5.1 and Eq.(LABEL:Eq._StochFloqBogolRelns)).

The TWA numerical calculations are based on expanding the stochastic field
functions in terms of time-independent gravitational modes, and during the
calculations $\widetilde{\psi}(z,t)$ is known at each time point (see Sect.
2.3.2 and Eqs. (48), (52), (53)). At each stage $\widetilde{\psi}(z,t)$ and
$\Phi_{c}(z,t)$ can also be re-expanded in terms of Floquet modes (see
Sect.2.3.1 and Eq.(45)).

A suitable description of the behaviour for both TWA and GPE numerics in terms
of the theoretical approach in this paper would be to think of a
representative stochastic field function or the condensate wave-function as
moving wave-packets, starting from their initial space-localised forms
$\widetilde{\psi}(z,0)$ and $\Phi_{c}(z,0)$ and then changing in shape and
position as time evolves. If $\widetilde{\psi}(z,t)$ or $\Phi_{c}(z,t)$ are
expanded in terms of Wannier modes, it may be the case that only one or two
modes are important, and the evolution of each component can be considered
separately. This situation can be described as two or more wave-packets that
may move in opposite directions, then combine again to reform the initial
$\widetilde{\psi}(z,t)$ or $\Phi_{c}(z,t)$ after a characteristic period has
elapsed.

As will be seen, there will be certain values for the boson-boson interaction
strength factor $gN$ that divide the observable behaviour between the presence
of a DTC and its absence, and also between where discrete time-translation
symmetry breaking occurs and where it does not. These changes occur in a
different way for Wannier and harmonic trap initial conditions; so to
distinguish the two cases, the values for $gN$ involved will be referred to as
the critical interaction strength or the threshold interaction strength,
respectively.

The numerical calculations are performed for a quasi-one-dimensional Bose-
Einstein condensate of ultracold atoms released from a harmonic trap to bounce
resonantly under the influence of gravity on an atom mirror oscillating with
period $T$ and amplitude $\lambda$ (in the oscillating frame). In carrying out
these calculations we use (dimensionless) gravitational units: length
$l_{0}=(\hbar^{2}/m^{2}g_{E})^{1/3}$, time
$t_{0}=(\hbar^{2}/mg_{E}^{2})^{1/3}$, energy $E_{0}=mg_{E}l_{0}$, where
$\hbar$ is the reduced Planck’s constant, $g_{E}$ is the Earth’s gravitational
acceleration and $m$ is the mass of the atom. For $s=2$ (period-doubling) we
set the drop height $h_{0}$ so that the bounce period
$t_{bounce}=(2h_{0})^{1/2}$ is equal to $2T$. In gravitational units the
potential in Eq. (3) is given by

$\overline{V}(\overline{z},\overline{t})=\overline{z}(1-\lambda\cos\overline{\omega}\overline{t})$
(104)

where $\overline{V}$ is the potential in units of $E_{0}$, and
$\overline{z},\overline{t}$ are position and time in units of $l_{0},t_{0}$
respectively, with $\overline{\omega}$ being the frequency in units of
$1/t_{0}$. Henceforth the over-bars will be left understood.

###  3.1 Choice of Parameters

As we wish to compare our TWA calculations with the mean-field GPE and time-
dependent Bogoliubov calculations in Kuros et al. [19], we chose the
parameters $\lambda$ and $\omega$ to be the same as theirs, namely
$\lambda=0.12$ and $\omega=1.4$ (in gravitational units). The choice made by
Kuros et al. [19] was informed by semi-classical treatments of the DTC system
[1]. These parameters allow a reasonably short atom transfer time (about $2$
$s$ for ${}^{7}Li$) between the two wave-packets for zero particle interaction
compared with a typical BEC lifetime, a reasonably small number of bounces
($<100$) during the atom transfer time, and an average interaction energy per
particle that is small compared with the energy gap between the first and
second bands in the single-band Bose-Hubbard model used in Ref. [19]. The
above values of $\lambda$ and $\omega$ result in two Floquet modes $\phi_{1}$
and $\phi_{2}$ at $t=0.5T$ that are similar to those in Fig. 2a of Ref. [1]
for $t=0T$. The parameters for the harmonic trap potential,
$\widetilde{\omega}_{0}=0.68$ and drop height, $h_{0}=9.82$, (in gravitational
units) for our TWA calculations were chosen so that at $t=0T$ the condensate
mode function has a large overlap with the Wannier function $\Phi_{2}$ that is
obtained from the two Floquet modes $\phi_{1}$ and $\phi_{2}$. The Floquet
frequencies are chosen to create a localised Wannier-like wave-packet which
leads to $(\nu_{2}-\nu_{1})$ being approximately $\omega/2$.

The various physical quantities for the simple case of a hard-wall mirror
potential $(V(z)=z$ for $z\geq 0$ and $\infty$ for $z<0)$ are summarised in
Table 1, both in gravitational units and in SI units for the case of the
bosonic ${}^{7}Li$ atom, which is chosen as an example because its $s$-wave
scattering length can be tuned precisely via a broad Feshbach resonance [18].
The drop height required to satisfy the resonance condition for $s=2$, i.e.,
$h_{0}l_{0}=20$ $\mu m$ for ${}^{7}Li$, is rather challenging to work with in
an experiment, but larger drop heights could be used by operating with larger
$s$ resonances [18]. The mirror oscillation amplitude, $\lambda
l_{0}/\omega^{2}=125$ $nm$, is for a theoretical hard-wall potential mirror;
in the case of a realistic soft Gaussian potential mirror the oscillation
amplitudes are typically an order of magnitude larger in order to achieve the
same driving effect [18].

Table 1

Physical quantities for $s=2$ resonance used in the calculations. All
quantities are in gravitational units except the expressions for $l_{0},t_{0}$
and $E_{0}.$

Quantity | Symbol | Gravitational | MKS units for ${}^{7}Li$
---|---|---|---
|  |  units |  (lab frame)
Gravitational unit of length | $l_{0}=(\hbar^{2}/m^{2}g_{E})^{1/3}$ |  |  $2.03$ $\mu$m
Gravitational unit of time | $t_{0}=(\hbar/mg_{E}^{2})^{1/3}$ |  |  $0.455$ ms
Gravitational unit of energy | $E_{0}=mg_{E}l_{0}$ |  |  $E_{0}/k_{B}=16.8$ nK
Mirror oscillation frequency | $\omega$ | $1.4$ |  $0.49$ kHz
Mirror oscillation period | $T=2\pi/\omega$ | $4.49$ |  $2.04$ ms
Mirror oscillation amplitude | $\lambda$ | $0.12$ |  $\lambda l_{0}/\omega^{2}=125$ nm
Bounce period for $s=2$ | $2T$ | $1.43$ |  $4.08$ ms
Drop height for $s=2$ | $h_{0}=\frac{1}{2}T^{2}$ | $9.82$ |  $20.0$ $\mu$m
Trap frequency | $\widetilde{\omega}_{0}$ | $0.68$ |  $238$ Hz
Boson transfer rate for $gN=0$ | $J$ |  $7.14\times 10^{-4}$ ∗ |  $1.6$ s-1
Boson transfer time for $gN=0$ | $t_{trans}=\pi/J$ | $440$ |  $2.0$ s
Number of bounces during $t_{trans}$ | $N_{b}=t_{trans}/T$ | $98$ | $98$
Longitudinal cloud width at $z=h_{0}$ | $\sigma_{z}=(\hbar/2m\omega_{0})^{1/2}/l_{0}$ | $0.85$ |  $1.7$ $\mu$m
s-wave scattering length | $a_{S}/a_{0}$ | $9.0\times 10^{3}|gN|/N$ | $9.0\times 10^{3}|gN|/N$

* Using mean -field calculations based on the overlap of two Wannier wave-functions used as a basis for the condensate wave-function, similar to calculations in Ref. [19].

###  3.2 Initial State - Linear Combination of Wannier-Like States

We first carry out TWA calculations for an initial condensate mode function
given by the $t=0$ form of the stable periodic solution of the time-dependent
mean-field GPE (Eq. (10)). This solution is approximated as a linear
combination of two single-particle Wannier-like modes $\Phi_{1}$ and
$\Phi_{2}$ (Eq. (59)), as in Ref. [1]. These Wannier-like modes each have
period $2T$ and are delayed by time $T$ with respect to each other, as shown
in Figs. 1 (c), (d), with the related Floquet modes shown in Figs. 1(a), (b).
The linear combination coefficients are obtained via the extrema of the energy
functional given by Eq. (3) in Ref. [1]. For$|gN|$ less than for the critical
interaction strength $-0.006$, the linear combination gives the Floquet mode
$\phi_{2}(z,t)$, which would have period $T$. When $|gN|\,>0.006$, a more
general linear combination of $\Phi_{1}(z,t)$ and $\Phi_{2}(z,t)$ results,
which would have period $2T$. Such a Wannier initial condition is useful for
studying the dynamical processes and demonstrating discrete time-translation
symmetry breaking, though it is not a realistic BEC state that can be prepared
in the laboratory. The position probability density (PPD) results are
presented in Fig. 2 and the one-body projector (OBP) and its Fourier transform
(FT) are shown in Fig. 3. The OBP for a many-body system, is similar to the
autocorrelation function or the fidelity used in Ref. [19] in a mean-field
calculation (see Sect. 2.2.8 for details).

The PPD plots in Fig. 2 are TWA (blue) and mean-field GPE (red) calculations
for interaction strengths $gN=-0.005$, $-0.006$ and $-0.007$, with $N=600$
bosonic atoms, and evolution times out to $t/T=1998$ mirror oscillations.
Initially, at time $t=0T$, the two wave-packets move towards each other; the
$\Phi_{1}$ wave-packet is reflected coherently from the oscillating mirror
creating fringes due to interference between the incident and reflected parts
of the wave-packet, while the $\Phi_{2}$ wave-packet is at the classical
turning point $z=h_{0}$. At times $t=(k+0.5)T$ $(k=0,1,2\ldots),$ the two
wave-packets cross paths and interfere, again creating fringes.

Figure 1: Single particle Floquet states and Wannier modes for $s=2$, initial
position $h_{0}=9.82$, trap frequency $\widetilde{\omega}_{0}=0.68$,
oscillation amplitude $\lambda=-0.12$ and drive frequency $\omega=1.4$. (a)
and (b) show the two Floquet states $\phi_{1}(z,t)$ and $\phi_{2}(z,t)$ with
the corresponding Floquet quasi-energies $\nu_{1}\approx 0.410$ and
$\nu_{2}\approx 1.109$. (c) and (d) show the two Wannier modes $\Phi_{1}(z,t)$
and $\Phi_{2}(z,t)$ constructed from the two Floquet states via Eq. (59)
Figure 2: Position probability density (PPD) TWA (blue) and mean-field GPE
(red) calculations and for Wannier initial conditions as a function of $z$ and
$t$ for $s=2$ and for different interaction strengths $gN=-0.005$, $-0.006$
and $-0.007$ in (a), (b) and (c), respectively. The calculations are carried
out for $N=600$, $h_{0}=9.82$, $\widetilde{\omega}_{0}=0.68$, $\lambda=0.12$
and $\omega=1.4$. Figure 3: One-body projector (OBP) and the corresponding
Fourier transform (FT) for Wannier initial conditions. (a), (b), (c) and (d)
show the OBP for $gN=-0.005$, $-0.006$, $-0.007$ and $-0.02$, respectively;
(e), (f), (g) and (h) shows the corresponding FT. The inset of (d) shows a
zoom-in of the OBP at a shorter time interval. Note that the pairs of OBP
versus $t$ plots are not continuous; the OBP oscillates between the two
curves, such as every $T$ as the inset for $gN=-0.02$ shows. The blue circle
(red cross) symbol show the TWA (GPE) results. Note the different $y$-axis
scale of (b) from (a), (c) and (d), as well as the different $y$-axis scale of
(e) and (f) from (g) and (h). The other parameters are the same as Figure 2.

For $gN=-0.005$ (Fig. 2 (a)), the PPDs are almost perfectly $T$-periodic and
the two Wannier wave-packets have about the same magnitude, while for
$gN=-0.007$ (Fig. 2 (c)), the PPDs are $2T$-periodic and the two wave-packets
have different magnitude. The change in response period of the BEC is clearly
illustrated in plots of the one-body projector versus evolution time $t/T$ and
its Fourier transform (Fig. 3). For $gN=-0.005,$ the OBP versus $t/T$ is
essentially a single oscillating horizontal line (Fig. 3 (a)), indicating a
period $T$ equal to the drive period, while for $gN=-0.006$, the OBP exhibits
a closely spaced double-curve structure (Fig. 3 (b)) and a small sub-harmonic
peak appears in the FT at half the driving frequency $\omega/2$ (Fig. 3 (f)).
Thus, discrete time-translation symmetry is broken at a critical interaction
strength near $gN\approx-0.006$ to form a DTC. For $gN=-0.007$ and $\
gN=-0.02$, the OBP exhibits a clear double-curve structure (Figs. 3(c), (d))
and a larger sub-harmonic peak in the FT at $\omega/2$ (Figs.3 (g), (h)),
indicating a single stable wave-packet for evolution times out to a least
$2000$ mirror oscillations.

For $gN=-0.006$ and $gN=-0.007$ (Figs. 2 (b), (c)), there is no noticeable
change in the PPDs between the set for $t=0$ to $2T$ and the set for $t=1996T$
to $1998T$, indicating the DTC shows no sign of decay for times out to at
least $t/T=1998$ mirror oscillations. Furthermore, there is no noticeable
difference between the TWA (red) and the mean-field GPE (blue) PPDs for both
$gN=-0.005$ and $gN=-0.007$. However, at the critical interaction strength,
$gN\approx-0.006$, there is a small difference in the OBP (up to $\approx
1\%$) between the TWA and the mean-field GPE calculations (Fig. 3 (b)) for
times longer than $500$ mirror oscillations, indicating a quantum depletion of
this order. This peak in the quantum depletion at the critical interaction
strength for discrete time-translation symmetry breaking is further
illustrated in Fig. 11.

In summary, in the case of the Wannier initial condition, for small $|gN|$ the
special linear combination of Wannier-like modes as the initial condition
evolves with a period $T$ as for a Floquet mode, while for $|gN|$ larger than
for the critical value of the interaction strength $gN\approx-0.006$, discrete
time-translation symmetry breaking occurs, and the more general linear
combination of Wannier-like modes now evolves with periodicity $2T$, to form a
DTC.

###  3.3  Initial State - Gaussian-like State Prepared in a Harmonic Trap

We now consider an initial Gaussian-like state that is prepared in a harmonic
trap and matches the Wannier wave-packet at the classical turning point. Such
an initial state can be realistically prepared in the laboratory. Our
calculations for the harmonic trap case are mainly for negative $gN.$

In Figs. 4-8, we present TWA (blue) and mean-field GPE (red) calculations of
the PPD for a range of interaction strengths $gN$, with $N=600$ bosonic atoms,
and evolution times out to $t/T=1998$ mirror oscillations. In Fig. 9 we show
the corresponding OBP and FT for the original condensate mode function and
Fig. 10 shows the OBP for the two Wannier modes $\Phi_{1}$, $\Phi_{2}$. Note
that the OBP for the original condensate mode involves a time-independent
function, whereas the OBP for the Wannier modes involve functions that are
time dependent.

For $gN=0$ (Fig. 4, based on the GPE), the PPD commences with $2T$-periodicity
- reflecting the bounce period - but after times $t=498T-500T$, about half the
bosonic atoms have transferred from the initial wave-packet to the second
wave-packet, while after times $t=998T-1000T$ essentially all of the atoms
have transferred to the second wave-packet, and then after times
$t=1996T-1998T$, all of the atoms have transferred back to the original wave-
packet, which is now a slightly asymmetric Wannier-like wave-packet (see Fig.
2). The transfer of bosonic atoms back and forth between the two wave-packets
appears as a periodic modulation of the OBP (Fig. 9 (a)) and a splitting of
the peak in the FT (Fig. 9 (e)). The magnitude of the splitting of the FT peak
matches a mean-field calculation of the coupling constant, $J=7.14\times
10^{-4}$, based on the overlap of the two Wannier wave-functions used as a
basis for the condensate wave-function, similar to calculations presented in
Ref. [19]. Confirmation that the bosonic atoms transfer back and forth between
the two wave-packets is provided by calculations of the number of atoms
occupying the first ($N_{1}$) and second wave-packets ($N_{2}$) versus time
(Fig. 10). Complete transfer of bosonic atoms between Wannier states
$\Phi_{1}$, $\Phi_{2}$ occurs at $gN=-0.006$ but the transfer starts to cease
at about $gN=-0.012$. Note that the PPD plot for $gN=-0.006$ (Fig. 5) is
almost identical to that for $gN=0$ (Fig. 4); so in Fig. 10 the $gN=-0.006$
plot also applies to $gN=0$. Also, in Fig. 5 for $gN=-0.006$ the TWA and GPE
plots of the PPD are indistinguishable, as they also would be for $gN=0$.

For $gN=-0.02$ (Fig. 7 - blue solid curve), the interaction is sufficiently
strong to suppress the transfer of bosonic atoms between the two wave-packets,
so that the initial wave-packet propagates as a single stable localised wave-
packet for times out to at least $t/T=1998$ mirror oscillations, indicating
the creation of a DTC. As may be seen in Fig. 10, this wave-packet is the
Wannier mode $\Phi_{2}$, which has a periodicity $2T$. The OBP (Fig. 9 (c))
further indicates little or no transfer of atoms between the two wave-packets
and the FT (Fig. 9 (g)) exhibits a single sub-harmonic peak at half the
driving frequency $\omega/2.$ Further evidence that there is little or no
transfer of atoms between the two wave-packets for $gN=-0.02$ is provided by
calculations of the atom numbers in the two wave-packets versus time (Fig.
10). For an even larger interaction strength $gN=-0.1$, the PPDs (Fig. 8) and
OBPs (Fig. 9 (d)) are similar to those for $gN=-0.02$ but now there is
essentially no transfer of atoms between the two wave-packets for times out to
at least $t/T=1998$ mirror oscillations. In addition, there is no significant
broadening of the wave-packet for times out to $t/T=1998$, indicating no
significant heating of the atom cloud by the periodic driving or by the
quantum many-body fluctuations. The creation of a single stable localised
wave-packet at such a large interaction strength indicates that the TWA
calculations are valid at interaction strengths up to at least $gN=-0.1$. This
is in a regime where for models based on a single-band description of the
Bose-Hubbard model [1], [17], [18] the interaction energy per particle becomes
comparable with the energy gap between the first and second bands and the
model is no longer valid. For $gN=0,-0.006,-0.02$ and $-0.1$, the TWA (blue)
and the mean-field GPE (red) PPDs are indistinguishable for times out to at
least $1998T$.

For $gN=-0.012$ (Fig. 6), the PPDs show only a weak transfer of bosonic atoms
between the two wave-packets, which is illustrated by the OBP and its FT
(Figs. 9 (b),(f)). Thus, $gN=-0.012$ represents the threshold interaction
strength for the creation of a single stable localised wave-packet (see Fig.
10). This represents the onset of DTC creation, since the atoms never
completely transfer out of the wave-packet $\Phi_{2}$, which is $2T$-periodic.
Furthermore, for times $t/T=1498$ $-1998$ mirror oscillations the wave-packet
for the TWA calculation is significantly smaller than for the mean-field GPE
calculation, indicating significant quantum depletion near this threshold
value. The difference (up to $\approx 30\%$) between the TWA and the mean-
field GPE calculations due to quantum depletion near the threshold interaction
strength is clearly illustrated in the OBP and its FT for times
$t/T=1498-1998$ mirror oscillations (Figs. 9 (b),(f)).

The red dashed PPD curves in Fig. 7 represent TWA calculations for a repulsive
interaction $gN=+0.02$. Interestingly, these are very similar to the PPD
curves for $gN=-0.02$ (blue solid curves), apart from a small shift in phase
of the interference fringes. In particular, there is no sign of decay or
broadening at times out to at least $1998$ mirror oscillations.

In summary, for harmonic trap initial conditions, we have a transition from a
wave-packet evolving with period $2T$ coupled (with coupling constant $J$ ) to
a second wave-packet with period $2T$, to a single stable wave-packet evolving
with period $2T$, and hence a DTC, at a threshold interaction strength
$gN\approx-0.012$. Below the threshold $gN$, the wave-packet has periodicity
$2T$, modulated at frequency $J$.

For the non-driving case with $gN=-0.10$ the PPD, OBP and its FT are presented
in Appendix 14 in the Supplementary Material (see Figs. 14, 15). No well-
defined periodicity is apparent in either the PPD or OBP results.

Figure 4: PPD for harmonic trap initial conditions as a function of $z$ and
$t$ for $s=2$ and non-interacting case ($gN=0$) for $h_{0}=9.82$,
$\widetilde{\omega}_{0}=0.68$, $N=600$, $\lambda=0.12$ and $\omega=1.4$ by
solving the time-dependent Schrödinger equation. Figure 5: PPD for harmonic
trap initial conditions as a function of $z$ and $t$ for interaction
$gN=-0.006$ using the same parameters as Figure 2. The blue solid (red dashed)
curves are calculated using the TWA (GPE) approach. Figure 6: The same as
Figure 5 but for $gN=-0.012$. The blue solid (red dashed) curves are
calculated using the TWA (GPE) approach. Figure 7: PPD for harmonic trap
initial conditions as a function of $z$ and $t$ for interaction $|gN|=0.02$
using the same parameters as Figure 2. The blue solid (red dashed) curves
correspond to the attractive (repulsive) interaction case. Figure 8: The same
as Figure 5 but for $gN=-0.1$. The blue solid (red dashed) curves are
calculated using the TWA (GPE) approach. Figure 9: .OBP and the corresponding
FT for harmonic trap initial conditions (a), (b), (c), and (d) show the OBP
for $gN=0$, $-0.012$, $-0.02$ and $-0.1$, respectively; (e), (f), (g) and (h)
show the corresponding FT. The red circle symbols in (a) and (e) show a
Schrödinger equation calculation for the non-interacting case. In other cases,
the blue circle (red cross) symbols show the TWA (GPE) results. The parameters
are the same as in Figure 2. Figure 10: Expectation value of particle numbers
$N_{1}$ and $N_{2}$ in single particle Wannier modes $\Phi_{1}$ and $\Phi_{2}$
as a function of time for different interaction strengths $gN$ shown in the
figure. The parameters are the same as in Figure 2.

###  3.4 Quantum Depletion

Figure 11 shows a logarithmic plot of the maximum quantum depletion due to
quantum many-body fluctuations obtained for the TWA calculations as a function
of both attractive and repulsive interaction strengths for a BEC of $N=600$
atoms and evolution times $t/T\leq 2000$ mirror oscillations. The calculations
were performed for both a Wannier-like initial state and a harmonic-trap
initial state with and without periodic driving of the mirror. Interestingly,
the plots are symmetrical about $gN=0$, as a result of the PPDs for an
attractive interaction being almost the same as those for a repulsive
interaction of the same magnitude (see Fig. 7).

For the case of a Wannier-like initial state, the quantum depletion in Fig. 11
is essentially zero for evolution times out to at least $t/T=2000$ mirror
oscillations, except close to the critical interaction strength for breaking
discrete time-translation symmetry, $|gN|=0.006$, where there are small peaks
in the quantum depletion corresponding to about $6$ atoms out of a total of
$N=600.$

For the case of a harmonic-trap initial state with driving, the quantum
depletion is less than $2$ atoms out of $N=600$ for evolution times out to at
least $t/T=2000$, except close to the threshold interaction strength for
creating a single localised wave-packet, $|gN|=0.012$, where the depletion is
as high as $260$ atoms out of $N=600$. To confirm the reliability of the TWA
calculations near the threshold interaction strength, we implemented an
analytical two-mode model for comparison and found excellent agreement within
the investigated time window. The two-mode model will be detailed elsewhere.
We note that our TWA calculations can handle a situation where the quantum
depletion is very large and where the time-dependent Bogoliubov theory
treatment of depletion [19] would break down.

For the case of a static mirror ($\lambda=0$), the TWA calculations indicate a
relatively large quantum depletion, i.e., $40-220$ atoms out of $N=600$ for
$|gN|\leq 0.1$. A comparison with the case of a resonantly driven mirror
indicates that the DTC created by resonant periodic driving together with a
sufficiently strong particle interaction significantly suppresses the quantum
depletion due to quantum many-body fluctuations.

Figure 12 (a) shows the quantum depletion and the occupation number of atoms
in the condensate mode ($p_{0}$) and the second occupied mode ($p_{1}$) versus
time out to $t/T=2000$ mirror oscillations for $gN=-0.015,N=600$, while Fig.
12 (b) shows the natural orbitals for the condensate mode and the second
occupied mode at$\;t/T=2000$. The total quantum depletion closely matches the
difference $N-p_{0}$ over the range $t/T=2000$, indicating that the depletion
is essentially due to the escape of atoms from the condensate mode 1 to the
second occupied mode 2. For $gN=-0.015$, the difference $N-p_{0}$ is very
close to the occupation number of the second occupied mode $p_{1}$, indicating
that the occupation number for all other modes is essentially zero ($<$
$1.6\%$) for times out to at least $t/T=2000$. This implies that models based
on just two modes should work well. For $gN=-0.015$, i.e., fairly close to the
threshold interaction strength for creating a single localised wave-packet,
the quantum depletion is about $5$ atoms out of $N=600$ for times out to
$t/T=1500$ and about $15$ atoms out of $N=600$ for times out to $t/T=2000$
mirror oscillations.

The quantum depletion versus $t/T$ curve in Fig. 12 suggests the quantum
depletion may be still increasing at $t/T=2000$. In Ref. [19] it is found that
in a two-mode theory treatment the quantum depletion for $gN=-0.02$ saturates
at times of about $t/T\approx 500\sqrt{N}$, or at $t/T\approx 12,000$ mirror
oscillations for $N=600$.

Figure 11: Maximum quantum depletion for $t<2000T$ as a function of $gN$. The
blue circles and purple squares correspond to the harmonic trap and Wannier
initial condition, respectively, where the calculation is carried out for
$N=600$ and other parameters are the same as Figure 2. We also show the
quantum depletion in red pentagrams for the harmonic trap initial condition
and a static mirror (no driving) for comparison. Figure 12: (a) Occupation
number as a function of time for harmonic trap initial condition with
$gN=-0.015$, $N=600$, and other parameters are the same as Figure 2. The blue
solid (orange dash-dotted) curve shows $N-p_{0}$ ($p_{1}$), respectively. We
also show the quantum depletion as red dashed curve. (b) The two mode
functions for the condensed mode (blue solid curve) and second occupied mode
(orange dashed curve) correspond to $p_{0}$ and $p_{1}$, respectively.

###  3.5 Mean Energy of Atoms at Long Times

Figure 13 presents TWA calculations of the mean energy $\left\langle
H\right\rangle$ for the case of a harmonic-trap initial state for times out to
$t/T=2000$ mirror oscillations and different interaction strengths $gN$. We
observe that the mean energy does not significantly increase and typically
oscillates around an average value close to its initial value. The evolution
of the mean energy indicates that the system reaches a steady state with no
net energy pumped from the drive and which is consistent with DTC behaviour.

Since our TWA method is a fully multi-mode approach that allows thermalisation
and occupations of many modes in the non-driving case with the same
interaction strength, we conclude that the absence of thermalisation in a DTC
is due to quantum effects of the resonant driving.

Figure 13: Mean energy as a function of time for harmonic trap initial
condition and different interaction strengths $gN$ shown in the figure. In
these calculations, $N=600$ and other parameters are the same as Figure 2.

##  4 Discussion

For a realistic initial condition corresponding to a harmonic trap condensate
mode function, our many-body TWA calculations for $s=2$ agree closely with the
mean-field GPE calculations for times out to at least $2000$ mirror
oscillations, except at interaction strengths $gN$ very close to the threshold
value for transition to a single stable wave-packet and DTC formation, where
the PPDs differ significantly from those determined from the mean-field theory