GPE approach. On the other hand, for a hypothetical initial condition in which
the BEC wavefunction is treated as a linear combination of Wannier-like
states, our TWA calculations agree broadly with the mean-field GPE
calculations, with only a small difference at interactions close to the
critical interaction for discrete time-translation symmetry breaking and DTC
formation.

For typical attractive interaction strengths greater in magnitude than the
threshold value for DTC formation and for the chosen parameters, the 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.
Here, our TWA calculations are also found to be in close agreement with recent
many-body calculations based on a time-dependent Bogoliubov approach [19].
Together with our results, this indicates that the occupation of additional
non-condensate modes is insignificant for most of the range of parameters
investigated. However, using the TWA approach we find that for interaction
strengths very close to the threshold value for DTC formation, the quantum
depletion due to the quantum fluctuations is as high as about $260$ atoms out
of a total of $600$ atoms, at times corresponding to $2000$ mirror
oscillations. Since Bogoliubov theory assumes that the non-condensate field is
relatively small, the calculations based on a time-dependent Bogoliubov
approach [19] would break down here, as would the PPDs calculations based on
mean-field theory [1], since the atoms are no longer in a single mode. The PPD
calculations based on the GPE differ from those obtained via the TWA approach,
as we have seen. For the range of conditions studied, the quantum depletion in
the vicinity of the threshold for creating a DTC results mainly from the
escape of atoms from the condensate mode to the second occupied mode and the
occupation of other modes is essentially zero ($<1.6\%$). This suggests that
many-body models based on just the two dominant modes should work well for
certain applications.

The general agreement we have found for our TWA calculations with those based
on the mean-field GPE and time-dependent Bogoliubov theory (TDBT) approaches
in Refs. [1], [17], [19], [18] (except near the threshold value for $gN$) is
not by chance. In Appendix 6 we show that the mean-field GPE and TDBT
equations can be derived from the phase-space TWA theory via an approximation
in situations where quantum depletion is small.

We also observe that the mean energy does not significantly increase for times
out to at least $2000$ mirror oscillations, 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 is consistent with a DTC behaviour. Thus, our TWA approach
predicts that thermalisation does not occur, at least out to $2000$ mirror
oscillations. Thermalisation cannot be treated using mean-field theory
approaches.

When the driving field is turned off and for a relatively strong interaction
($0.05\leq|gN|\leq 0.1$), the interaction is strong enough to couple many
modes and we find using the TWA that the quantum depletion can be more than
$220$ atoms out of a total of $600$ atoms at times corresponding to $2000$
mirror oscillations. However, with driving, the quantum depletion associated
with quantum fluctuations is strongly suppressed, and thermalisation quenched,
implying that the absence of our system’s thermalisation is a genuine many-
body effect, and along with DTC formation, is a direct consequence of driving
in the presence of a sufficiently strong interaction.

Thus, in certain circumstances, namely for realistic harmonic trap initial
conditions and where the interaction parameter $gN$ is close to the threshold
value for DTC formation ($gN=-0.012$ for our parameters), previous theories
based on mean-field theory ([1], [19]) or time-dependent Bogoliubov theory
[19] break down, and the phase-space theory TWA approach is required in order
to provide a comprehensive treatment of DTC behaviour covering all regimes.

We also find that the dynamical behaviour of our system is largely independent
of whether the boson-boson interaction is attractive or repulsive, and that
for sufficiently large repulsive interactions it is possible to create a
stable DTC based on repulsive interactions. The ability to create a stable
time crystal with repulsive interactions should allow more flexibility; for
example, the use of repulsive interactions may allow us to create a DTC at
still larger interaction strengths where there is no limitation imposed by
bosenova collapse [40] or the formation of a bright soliton as there is with
attractive interactions.

The explanation of the behaviour being essentially the same for negative and
positive $g$ can be provided from the Ito stochastic field equations (8), (9),
noting that all the observable quantities are determined from
$\widetilde{\psi}(z,t)$ and $\widetilde{\psi}^{+}(z,t)$. For large $N$ the
$\delta_{C}(z,z)$ term can be ignored, and the solutions for $-g$ can be
obtained from those for $+g$ via the substitution
$\widetilde{\psi}(z,t)=i\,\widetilde{\psi}(z,t)^{\\#},\widetilde{\psi}^{+}(z,t)=i\,\widetilde{\psi}^{+}(z,t)^{\\#}$
\- which is possible for the double phase-space Wigner distribution $W^{+}$
where the field functions and stochastic field functions
$\widetilde{\psi}(z,t)$ and $\widetilde{\psi}^{+}(z,t)$ are allowed to be
unrelated. The PPD, etc., calculated from the new fields would be for $-g$,
and would have the same $z,t$ behaviour as for $+g$.

##  5 Summary and Conclusions

We have presented a full theoretical many-body quantum study based on
physically realisable initial conditions for creating discrete time crystals
using a Bose-Einstein condensate which is allowed to bounce resonantly on an
oscillating mirror in a gravitational field, as proposed in [1]. Our theory of
DTC creation allows for the effects of there being many modes that the bosons
could occupy; it allows for the mutually important effects of both
interactions and driving; it determines whether or not thermalisation occurs;
and it enables the sensitivity of DTC behaviour to changes in the initial
conditions to be studied - including the effect of the initial BEC
temperature.

The significance of the theory presented here is that it is not based on
restricted assumptions that have applied to previous work on DTC formation in
such systems, such as a mean-field theory ([1] and [19]), time-dependent
Bogoliubov theory ([19]) or a two-mode model. In a mean-field theory, all the
bosons are assumed to remain in a single mode, and a time-dependent Bogoliubov
theory is restricted to situations where the depletion from the condensate
mode is small. In contrast, our TWA theory allows for the large depletion and
quantum fluctuations that occur near the threshold $gN$ value for a DTC to
form. Thermalisation can also be treated, unlike in mean-field theory.

The ability to create a DTC over a broad range of attractive and repulsive
interaction strengths also provides flexibility for applications of DTCs; for
example, to the study of a wide range of non-trivial condensed matter
phenomena in the time domain. Condensed matter phenomena that have been
proposed based on a driven bouncing BEC platform include Anderson localization
[30] and many-body localization [35] in the time domain due to temporal
disorder; Mott insulator-like phases in the time domain [30]; quasi-
crystalline structures in time [36]; topological time crystals [37]; time
crystals with exotic long-range interactions [38]; and dynamical quantum phase
transitions in time crystals [39].

Possible future directions of this work include applying the TWA approach for
times much longer than $2000T$. We propose to investigate whether long-time
calculations could be based on using relatively small numbers of Wannier
modes, rather than the large number of gravitational modes used in the present
study. We also propose to apply the TWA approach to treat higher order
$sT$-periodicity, such as in the range $s=5-10$. This type of DTC has been
predicted using mean-field theory ([17], [18]), but the conditions for
creating a DTC could be investigated using a theoretical approach that does
not rely on the mean-field assumption. The TWA theory is already formulated
for treating finite temperature effects, and could be used to investigate the
conditions under which a DTC can be created at non-zero temperatures. Finally,
we wish to explore more fully the conditions for $sT$-periodicity than is
possible in this initial paper, such as the range of the boson-boson
interaction strength, different initial conditions and the stability of the
DTC. The upper bound of useable attractive interaction strengths is limited by
bosenova collapse [40], and it would be interesting to also study this regime.

###  5.1 Acknowledgements

We thank Krzysztof Sacha for many fruitful discussions and suggestions, and
for a critical reading of this article. Also, Peter Drummond, King Lun Ng and
Run Yan Teh are acknowledged for discussions on phase-space theory. Support of
the Australian Research Council (DP190100815) is gratefully acknowledged. JW
acknowledges an ARC Discovery Early Career Researcher Award (DE180100592).

###  5.2 Orcid ID

These are:

Jia Wang 0000-0002-9064-5245

Peter Hannaford 0000-0001-9896-7284

Bryan Dalton 0000-0002-1176-7528

##  6 Appendix A - Mean-Field Theory and Time-Dependent Bogoliubov Theory

In this appendix we show how the mean-field Floquet theory equation set out as
Eq. (2) in Ref [1], and which is the key starting equation in subsequent
papers [17], [19], [18], can be derived as an approximation from the phase-
space truncated Wigner approximation (TWA) theory. As we will see, Eq (2) is
the equation that a condensate field function contribution in the stochastic
field function would satisfy. We also show how the time-dependent Bogoiliubov
theory, which is a further development in Ref. [19] can also be seen as an
approximation to the TWA approach. The starting point is the Ito stochastic
field equation for the stochastic field function, but now based on the grand
canonical Hamiltonian $\widehat{H}_{G}=\widehat{H}-\mu\widehat{N}$, where
$\widehat{H}$ is the Hamiltonian, $\widehat{N}$ is the number operator and
$\mu$ is the chemical potential. We only consider number conserving states. As
usual we consider a regime where the boson number $N$ is large, $N\gg 1$. The
advantage of this approach is that Eq (2) is now based on a standard many-body
theory approach, whereas at present its justification rests on it being a
physically reasonable combination of the Gross-Pitaeskii equation (GPE) for a
single-mode condensate and the Floquet equation for a single-particle system
in a periodic potential.

###  6.1 FFPE and Ito SFE - Extra Term

The density operator $\widehat{\rho}$ satisfies the Liouville-von Neumann
equation

$\frac{\partial}{\partial
t}\widehat{\rho}=\frac{1}{i\hbar}[\widehat{H},\widehat{\rho}]$ (105)

but for number conserving states $[\widehat{N},\widehat{\rho}]=0$, so that we
also have

$\frac{\partial}{\partial
t}\widehat{\rho}=\frac{1}{i\hbar}[\widehat{H}_{G},\widehat{\rho}]$ (106)

involving the grand canonical Hamitonian.

The functional Fokker-Planck equation (see Eqs. (4),
(LABEL:Eq.FFPEPotentialFermiFldModel) and (6)) for the Wigner distribution
functional $W$ then has an additional term

$\left(\frac{\partial}{\partial
t}W\right)_{\mu}=-\frac{\mu}{i\hbar}\left(-\frac{\delta}{\delta\psi}(\psi
W)+\frac{\delta}{\delta\psi^{+}}(\psi^{+}W)\right)$ (107)

and this results in an extra term in the Ito stochastic field equations on the
right hand side

$\left(\frac{\partial\widetilde{\psi}}{\partial
t}\right)_{\mu}=\frac{1}{i\hbar}(-\mu\,\widetilde{\psi})\qquad\qquad\left(\frac{\partial\widetilde{\psi}^{+}}{\partial
t}\right)_{\mu}=\frac{1}{i\hbar}(+\mu\,\widetilde{\psi}^{+})$ (108)

Consequently, from Eqs. (8) and (9), the full Ito stochastic field equations
are now

$\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)-\mu\widetilde{\psi}(z,t)\right]$

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)-\mu\widetilde{\psi}^{+}(z,t)\right]$

###  6.2 Condensate and Non-Condensate Components

We first make the large $N$ approximation of neglecting the terms
$\delta_{C}(z,z)$ in comparison to
$\widetilde{\psi}^{+}(z,t)\widetilde{\psi}(z,t)$ since the former term is of
order $1$ whereas the latter is of order $N$.

We then write the stochastic field functions as the sum of a deterministic
term $\Phi_{0}(z,t)$ (or its conjugate) corresponding to the condensate mode
and a stochastic term $\delta\widetilde{\psi}(z,t)$ (or
$\delta\widetilde{\psi}^{+}(z,t)$) corresponding to the non-condensate modes.
The stochastic terms $\delta\widetilde{\psi}(z,t)$ or
$\delta\widetilde{\psi}^{+}(z,t)$ are assumed to be small. A similar approach
is made in time-dependent Bogoliubov theory in regard to the field operators.

Thus, we have

$\widetilde{\psi}(z,t)=\Phi_{0}(z,t)+\delta\widetilde{\psi}(z,t)\qquad\widetilde{\psi}^{+}(z,t)=\Phi_{0}^{\ast}(z,t)+\delta\widetilde{\psi}^{+}(z,t)$
(111)

If we substitute into Eqs. (LABEL:Eq.ItoSFEAMu) and (LABEL:Eq.ItoSFABMu) and
retain only the lowest order terms we find that

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

This is the same as Eq.(2) in Ref. [1].

The stochastic non-condensate terms satisfy the equations

$\displaystyle\frac{\partial}{\partial t}\delta\widetilde{\psi}(z,t)=$
$\displaystyle-\frac{i}{\hslash}\left[-\frac{\hslash^{2}}{2m}\frac{\partial^{2}}{\partial
z^{2}}+V(z,t)-\mu+2g\,\Phi_{0}^{\ast}(z,t)\Phi_{0}(z,t)\right]\;\delta\widetilde{\psi}(z,t)$
$\displaystyle-\frac{i}{\hslash}\left[g\,\Phi_{0}(z,t)^{2}\right]\;\delta\widetilde{\psi}^{+}(z,t)$
(113)

and

$\displaystyle\frac{\partial}{\partial t}\delta\widetilde{\psi}^{+}(z,t)=$
$\displaystyle+\frac{i}{\hslash}\left[-\frac{\hslash^{2}}{2m}\frac{\partial^{2}}{\partial
z^{2}}+V(z,t)-\mu+2g\,\Phi_{0}^{\ast}(z,t)\Phi_{0}(z,t)\right]\;\delta\widetilde{\psi}^{+}(z,t)$
$\displaystyle+\frac{i}{\hslash}\left[g\,\Phi_{0}^{\ast}(z,t)^{2}\right]\;\delta\widetilde{\psi}(z,t)$
(114)

correct to the lowest order in $\delta\widetilde{\psi}(z,t)$ and
$\delta\widetilde{\psi}^{+}(z,t)$. Here, we see that the stochastic
fluctuation fields $\delta\widetilde{\psi}(z,t)$ and
$\delta\widetilde{\psi}^{+}(z,t)$ are coupled. The last two equations are
analogous to time-dependent Bogoliubov-de Gennes equations (compare the time-
independent Eqs. (77)). Note the well-established $2g$ factor that appears in
Bogoliubov theory. The difference is that within the TWA approach the non-
condensate fields $\delta\widetilde{\psi}(z,t)$ and
$\delta\widetilde{\psi}^{+}(z,t)$ are treated as stochastic quantities,
whereas in time-dependent Bogoliubov theory, the equivalent quantities would
be time-dependent Bogoliubov mode functions $\overline{u}_{k}(z,t)$,
$\overline{v}_{k}(z,t)$.

The conclusion then is that the approach to mean-field theory and time-
dependent Bogoliubov theory in Refs. [1], [17], [19], [18] are now shown to be
an approximation to the TWA theory, valid in regimes where almost all the
bosons are in one condensate mode and the quantum depletion into non-
condensate modes is relatively small. There will of course be regimes where
this approximate solution to the TWA equations will break down, and in future
work we hope to explore some of these.

###  6.3 Energy Functional

In Ref. [1] Eq.(2) can be derived by minimising the energy functional set out
in Eq.(3) of Ref. [1], subject to the normalisation constraint that
${\displaystyle\int}dz\;\Phi_{0}^{\ast}(z,t)\Phi_{0}(z,t)=N.$ The energy
functional involves an integral over time $t$, which in Ref. [1] is taken over
the interval $0$ to $2T$, where $T$ is the period of the driving potential.
The energy functional could be taken as an average over interval $sT$ in
general

$E[\Phi_{0}]=\frac{1}{sT}{\displaystyle\int}dz{\displaystyle\int\limits_{0}^{sT}}dt\;\Phi_{0}^{\ast}\left[-\frac{\hslash^{2}}{2m}\frac{\partial^{2}}{\partial
z^{2}}+V(z,t)-i\hbar\frac{\partial}{\partial
t}+\frac{g}{2}\Phi_{0}^{\ast}\Phi_{0}\right]\Phi_{0}$ (115)

so minimising this with respect to changes in $\Phi_{0}$ subject to the
contraint leads to Eq. (LABEL:Eq.GPEFloquet) above for $\Phi_{0}$. The
chemical potential appears via a Lagrange undetermined multiplier. The choice
of the time interval for a different $s$ does not result in a different
equation for $\Phi_{0}$.

The advantage of this linkage is that if the condensate wave-function is
expanded in terms of Wannier functions $\Phi_{i}(z,t)$ which have periodicity
$sT$

$\Phi_{0}(z,t)={\displaystyle\sum\limits_{i}}a_{i}(t)\Phi_{i}(z,t)$ (116)

the energy functional becomes a function $E(a_{i})$ of the coefficients
involving a fourth-order polynomial. This can then be minimised to give the
condensate wave-function. The energy function will then involve coefficients
that are space-time integrals involving the Wannier functions up to the fourth
power. However, there are no extra exponential factors of the form
$\exp(-in\omega t)$ that occur in many-body formulations of Floquet theory
based on Shirley’s paper [41], and which can be eliminated via a fast-rotating
approximation.

##  7 Appendix B - Details regarding TWA Validity

Some specific questions regarding the TWA validity and reliability are as
follows:

(a) Does the accuracy of the TWA calculations change when the interactions
change from weak to strong? Does its accuracy depend on the extent of
correlations between different modes of the system?

At present we have no reason to believe that our TWA calculations would not be
accurate for the weak interaction situation we studied. Even if (as according
to Ref. [21]) the TWA is not suitable to describe systems with very strong
interaction and strong correlations, in our system we focus on the weak
interaction situation. We strongly believe though that the TWA can treat the
strong interaction regime, as the derivation of the final Ito stochastic field
equations does not depend on the size of the coupling constant $g$. In regard
to correlations, our multi-mode TWA approach allows for correlations between
modes to develop, so the TWA should be able to treat strong correlations. Many
gravitational modes become correlated during the evolution, but for our weak
interaction case the main correlation is between only two Wannier modes.

(b) Do the TWA predictions hold for long evolution times or will build up of
errors become uncontrollable? Is the time over which the TWA remains accurate
long enough to confirm time crystal behaviour?

The evolution error is a numerics issue rather than one involving the theory’s
formalism. In principle, the TWA method can accumulate error in our adopted
Runge-Kutta numerical method as the evolution time increases. However, in our
chosen time-window $(t<2000T)$, convergent tests with respect to time steps
have been performed and show very good convergence. Stable time-crystal
behaviour has been seen within this time window for $|gN|>0.012$ with the
harmonic trap initial condition. The time evolution for even longer times is
currently limited by computational resources. As indicated in the last
paragraph of the Summary, we plan a further investigation using more advanced
computers regarding whether the time crystal will remain stable after much
longer evolution times.

(c) What will the effect of the TWA approximation be on the results? Which
features will display differences with respect to exact behaviour?

Estimating the error coming from the approximation that neglects the third-
order functional derivative term is very difficult. As far as we know, no
exact calculation exists for the solution of the full FFPE for the sort of
many-boson systems we are studying and no equivalent Ito stochastic field
equations have even been derived when third-order derivative terms are
included. Therefore, this important question will remain unanswered until
researchers find a way to obtain exact solutions for the full FFPE. This is
beyond the scope of the present paper.

## References

  * [1] K. Sacha. Modeling spontaneous breaking of time-translation symmetry. Phys. Rev. A 91, 033617 (2015).
  * [2] V. Khemani, A. Lazarides, R. Moessner and S. L. Sondhi. Phase structure of driven quantum systems. Phys. Rev. Lett. 116, 250401 (2016).
  * [3] D. V. Else, B. Bauer and C. Nayak. Floquet time crystals. Phys. Rev. Lett. 117, 090402 (2016).
  * [4] N. Y. Yao, A. C. Potter, I. D. Potirniche and A. Vishwanath. Discrete time crystals. rigidity, criticality, and realizations. Phys. Rev. Lett. 118, 030401 (2017).
  * [5] K. Sacha. Time Crystals (Berlin: Springer, 2020).
  * [6] J. Zhang, P. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I. D. Potirniche, A. C. Potter, A. Vishwanath et al.. Observation of a discrete time crystal. Nature 543, 217 (2017).
  * [7] S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani et al.. Observation of dicrete time-crystalline order in a disordered dipolar many body system. Nature 543, 221 (2017).
  * [8] S. Pal, N. Nishad, T. S. Mahesh and G. J. Sreejith. Temporal order in periodically driven spins in star-shaped clusters. Phys. Rev. Lett. 120, 180602 (2018).
  * [9] J. Rovny, R L Blum and S E Barrett. Observation of discrete-time-crystal signatures in an ordered dipolar many-body system. Phys. Rev. Lett. 120, 180603 (2018).
  * [10] J. Rovny, R. L. Blum and S. E. Barrett. ${}^{31}P$ NMR study of discrete time-crystalline signatures in an ordered crystal of ammonium dihydrogen phosphate. Phys. Rev. B 97, 184301 (2018).
  * [11] J. Smits, L. Liao, H. T. C. Stoof and P. van der Straten. Observation of a space-time crystal in a superfluid quantum gas. Phys. Rev. Lett. 121, 185301 (2018).
  * [12] L. Liao, J. Smits, P. van der Straten and H. T. C. Stoof. Dynamics of a space-time crystal in an atomic Bose-Einsten condensate. Phys. Rev. A 99, 013625 (2018).
  * [13] J. Smits, H. T. C. Stoof and P. van der Straten. On the long-term stability of space-time crystals. ArXiv:2007.07038 (2020).
  * [14] K. Sacha and J. Zakrzewski. Time crystals. Rep. Prog. Phys. 81, 016401 (2018).
  * [15] N. Y. Yao and C. Nayak. Time crystals in periodically driven systems. Phys. Today Sept Issue 44 (2018).
  * [16] V. Khemani, R. Moessner and S. L. Sondhi. A brief history of time crystals. ArXiv 1910.10745 (2019).
  * [17] K. Giergiel, A. Kosior, P. Hannaford and K. Sacha. Time crystals: analysis of experimental conditions. Phys. Rev. A 98, 013613 (2018).
  * [18] K. Giergiel, T. Tran, A. Zaheer, A. Singh, A. Sidorov, K. Sacha and P. Hannaford. Creating big time crystals with ultracold atoms, New J. Phys. 22, 085004 (2020).
  * [19] A. Kuroś, R. Mukherjee, W. Golletz, F. Sauvage, K. Giergiel, F. Mintert and K. Sacha. Phase diagram and optimal control for n-tupling discrete time crystal. New J. Phys. 22, 095001. (2020).
  * [20] M. J. Steel, M. K. Olsen, L. I. Plinek, P. D. Drummond, S. M. Tan, M. J. Collet, D. F. Walls and R. Graham. Dynamical quantum noise in trapped Bose-Einstein condensates, Phys. Rev. A 58, 4824 (1998).
  * [21] P. B. Blakie, A. S. Bradley, M. J. Davis, R. J. Ballagh and C. W. Gardiner. Dynamics and statistical mechanics of ultra-cold Bose gases using c-field techniques, Adv. Phys. 57, 363 (2008).
  * [22] B. J. Dalton, J. Jeffers and S. M. Barnett. Phase Space Methods for Degenerate Quantum Gases (Oxford: Oxford University Press, 2015).
  * [23] C. W. Gardiner and P. Zoller. Quantum World of Ultra-Cold Atoms and Light: Book 3 - Ultra-Cold Atoms (Singapore: World Scientific, 2017).
  * [24] K. L. Ng, R. Polkinghorne, B. Opanchuk, and P. D. Drummond. Phase-space representations of thermal Bose-Einstein condensates, J. Phys. A: Math. Theor. 52, 035302 (2019).
  * [25] S. A. Morgan. A gapless theory of Bose-Einstein condensation in dilute gases at finite temperature, J. Phys B: At. Mol. Opt. Phys. 33, 3847 (2000).
  * [26] M. Lewenstein and L. You. Quantum phase diffusion of a Bose-Einstein condensate, Phys. Rev. Lett. 77, 3489 (1996); P. Villain, M. Lewenstein, R. Dum, Y. Casti, L. You, A. Imamoglu and T. A. B. Kennedy. Quantum dynamics of the phase of a Bose-Einstein condensate. J. Mod. Opt. 44, 1775 (1997).
  * [27] N. P. Proukakis and B. Jackson. Finite-temperature models of Bose-Einstein condensation. J. Phys. B: At. Mol. Opt. Phys. 41, 203002 (2008).
  * [28] A. J. Leggett. Bose-Einstein condensation in the alkali gases: Some fundamental concepts. Rev. Mod. Phys. 73, 307 (2001).
  * [29] J. H. Shirley. Solution of the Schrödinger equation with a Hamitonian periodic in time. Phys. Rev. 138, B979 (1965).
  * [30] K. Sacha Anderson localisation and Mott insulator phase in the time domain. Sci. Rep. 5 10787 (2015).
  * [31] N. Bogoliubov. On the theory of superfluidity. J. Phys. (USSR) 11, 23 (1947).
  * [32] A. Fetter. Non uniform states of an imperfect Bose gas. Ann. Phys. 70 67 (1972).
  * [33] G. C. Wick. The evolution of the collision matrix. Phys. Rev. 80, 268 (1950).
  * [34] M. K. Olsen, A. S. Bradley and S. B. Cavalcanti. Fock-state dynamics in Raman photoassociation of Bose-Einstein condensates. Phys Rev A 70, 033611 (2004).
  * [35] M. Mierzejewski, K. Giergiel and K. Sacha. Many-body localization caused by temporal disorder. Phys. Rev. B 96, 140201 (2017).
  * [36] K. Giergiel, A. Kuros and K. Sacha. Discrete time quasicrystals. Phys. Rev. B 99, 220303 (2019).
  * [37] K. Giergiel, A. Dauphin, M. Lewenstein, J. Zakrzewski and K. Sacha. Topological time crystals, New J. Phys. 21, 052003 (2019).
  * [38] K. Giergiel, A. Miroszewski and K. Sacha. Time crystal platform: From quasicrystal structures in time to systems with exotic interactions, Phys. Rev. Lett. 120, 140401 (2018).
  * [39] A. Kosior and K. Sacha. Dynamical quantum phase transitions in discrete time crystals. Phys. Rev. A 97, 053621 (2018).
  * [40] J. L. Roberts, N. R. Claussen, S. L. Cornish, E. A. Donley, E. A. Cornell and C. E. Wieman. Controlled collapse of a Bose-Einstein condensate, Phys. Rev. Lett. 86, 4211 (2001).
  * [41] J Shirley. Solution of the Schrodinger Equation with a Hamiltonian Periodic in Time. Phys. Rev. 165, B480 (1965).

SUPPLEMENTARY
MATERIAL

##  8 Appendix S1 - Mean Energy

###  8.1 General Expression

To derive the TWA stochastic expression (13) for the mean energy Eq. (12) for
the Hamiltonian in Eq. (1) we must first replace the _normally ordered_
products involving the field operators $\hat{\Psi}(z)$,
$\hat{\Psi}(z)^{{\dagger}}$ in each of the kinetic energy $\widehat{K}$,
potential energy $\widehat{V}$ and interaction energy $\widehat{U}$ terms by
their _symmetrically ordered_ forms (see for example Sect 7.1.2 in Ref [22]).
This can be done using Wick’s theorem [33]. Expansions of the field operators
in terms of a suitable set of orthonormal mode functions $\phi_{k}(z)$ are
also used.

We find that in the kinetic energy term

$\displaystyle\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)$ $\displaystyle=\left\\{\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)\right\\}-\frac{1}{2}\triangle\delta_{C}(z,z)$ (117)

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}$ (118)

In the potential energy term we have

$\hat{\Psi}(z)^{{\dagger}}\hat{\Psi}(z)=\left\\{\hat{\Psi}(z)^{{\dagger}}\hat{\Psi}(z)\right\\}-\frac{1}{2}\delta_{C}(z,z)$
(119)

The interaction energy term is more difficult to deal with, but in this term
we find that

$\displaystyle\hat{\Psi}(z)^{{\dagger}}\hat{\Psi}(z)^{{\dagger}}\hat{\Psi}(z)\hat{\Psi}(z)$
$\displaystyle=\left\\{\hat{\Psi}(z)^{{\dagger}}\hat{\Psi}(z)^{{\dagger}}\hat{\Psi}(z)\hat{\Psi}(z)\right\\}$
$\displaystyle-2\,\delta_{C}(z,z)\,\left\\{\hat{\Psi}(z)^{{\dagger}}\hat{\Psi}(z)\right\\}$
$\displaystyle+\frac{1}{2}(\delta_{C}(z,z))^{2}$ (120)

Hence, overall we have for the Hamiltonian

$\displaystyle\widehat{K}$
$\displaystyle=\int\mbox{\rule{-2.84526pt}{0.0pt}}dz\,\frac{\hbar^{2}}{2m}\left(\left\\{\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)\right\\}-\frac{1}{2}\triangle\delta_{C}(z,z)\right)$
$\displaystyle\widehat{V}$
$\displaystyle=\int\mbox{\rule{-2.84526pt}{0.0pt}}dz\,V(z,t)\left(\left\\{\hat{\Psi}(z)^{{\dagger}}\hat{\Psi}(z\right\\}-\frac{1}{2}\delta_{C}(z,z)\right)$
$\displaystyle\widehat{U}$
$\displaystyle=\int\mbox{\rule{-2.84526pt}{0.0pt}}dz\,\frac{g}{2}\left(\left\\{\hat{\Psi}(z)^{{\dagger}}\hat{\Psi}(z)^{{\dagger}}\hat{\Psi}(z)\hat{\Psi}(z)\right\\}-2\,\delta_{C}(z,z)\,\left\\{\hat{\Psi}(z)^{{\dagger}}\hat{\Psi}(z)\right\\}+\frac{1}{2}(\delta_{C}(z,z))^{2}\right)$
(121)

The final step is to use the results for the Wigner distribution functional
giving the mean value for a symmetrically ordered forms involving the field
operators

$\displaystyle\left\langle\left\\{\hat{\Psi}(z^{\\#})^{{\dagger}}\hat{\Psi}(z)\right\\}\right\rangle$
$\displaystyle={\displaystyle\int}D^{2}\psi\,D^{2}\psi^{+}\;\psi(z)\,\psi^{+}(z^{\\#})\,W[\psi,\psi^{+}]$
$\displaystyle=\overline{\widetilde{\psi}(z,t)\,\widetilde{\psi}^{+}(z^{\\#},t)}\,$
$\displaystyle\left\langle\left\\{\partial_{z}\hat{\Psi}(z^{\\#})^{{\dagger}}\,\partial_{z}\hat{\Psi}(z)\right\\}\right\rangle$
$\displaystyle={\displaystyle\int}D^{2}\psi\,D^{2}\psi^{+}\;\partial_{z}\psi(z)\,\partial_{z}\psi^{+}(z^{\\#})\,W[\psi,\psi^{+}]$
$\displaystyle=\overline{\partial_{z}\widetilde{\psi}(z,t)\,\partial_{z}\widetilde{\psi}^{+}(z^{\\#},t)}$
$\displaystyle\left\langle\left\\{\hat{\Psi}(z_{1}^{\\#})^{{\dagger}}\hat{\Psi}(z_{2}^{\\#})^{{\dagger}}\hat{\Psi}(z_{1})\hat{\Psi}(z_{2})\right\\}\right\rangle$
$\displaystyle={\displaystyle\int}D^{2}\psi\,D^{2}\psi^{+}\;\psi(z_{1})\,\psi(z_{2})\,\psi^{+}(z_{1}^{\\#})\,\psi^{+}(z_{2}^{\\#})\,W[\psi,\psi^{+}]$
$\displaystyle=\overline{\widetilde{\psi}(z_{1},t)\,\widetilde{\psi}(z_{2},t)\,\widetilde{\psi}^{+}(z_{1}^{\\#},t)\,\widetilde{\psi}^{+}(z_{2}^{\\#},t)}$
(122)

where we have expressed the mean value of symmetrically ordered forms
involving the field operators first as phase space functional integrals
involvinng the field functions $\psi(z),\psi^{+}(z)$ then as stochastic
averages involving the stochastic field functions
$\widetilde{\psi}(z,t),\widetilde{\psi}^{+}(z,t)$ (see Sects 12.3.3 and 15.1.9
in Ref [22]). The final stochastic average expression for
$\left\langle\widehat{H}\right\rangle$ is given in Eq (13).

###  8.2 Gravitational Modes

The mean energy expression (13) can be simplified via expanding the stochastic
field functions $\widetilde{\psi}(z,t),\widetilde{\psi}^{+}(z,t)$ in some
terms via gravitational modes.

The potential energy is $V(z,t)=mg_{E}z-mg_{E}\lambda z\cos\omega t$ so the
mean value of $\widehat{T}+\widehat{V}$ is

$\displaystyle\left\langle\widehat{T}+\widehat{V}\right\rangle$
$\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+{\displaystyle\int}dz\,\overline{\widetilde{\psi}^{+}(z,t)\left(-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial
z^{2}}+mg_{E}z\right)\,\widetilde{\psi}(z,t)}$
$\displaystyle-\frac{1}{2}{\displaystyle\int}dz\,\left(\frac{\hbar^{2}}{2m}\Delta\delta_{C}(z,z)-mg_{E}z\delta_{C}(z,z)\right)$
(123)

In the second and third lines we expand
$\widetilde{\psi}(z,t),\widetilde{\psi}^{+}(z,t)$ using Eq.(48) and
$\Delta\delta_{C}(z,z),\delta_{C}(z,z)$ via (7) and (14) using gravitational
modes. This gives after spatial integration by parts

$\displaystyle{\displaystyle\int}dz\,\overline{\widetilde{\psi}^{+}(z,t)\left(-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial
z^{2}}+mg_{E}z\right)\,\widetilde{\psi}(z,t)}$
$\displaystyle-\frac{1}{2}{\displaystyle\int}dz\,\left(\frac{\hbar^{2}}{2m}\Delta\delta_{C}(z,z)-mg_{E}z\delta_{C}(z,z)\right)$
$\displaystyle={\displaystyle\sum\limits_{k,l}}\overline{\widetilde{\eta}_{k}^{+}(t)\widetilde{\eta}_{l}(t)}\,{\displaystyle\int}dz\,\xi_{k}^{\ast}(z)\,\left(-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial
z^{2}}+mg_{E}z\right)\xi_{l}(z)$
$\displaystyle-\frac{1}{2}{\displaystyle\int}dz\,{\displaystyle\sum\limits_{k}}\xi_{k}^{\ast}(z)\,\left(-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial
z^{2}}-mg_{E}z\right)\xi_{k}(z)$
$\displaystyle={\displaystyle\sum\limits_{k,}}\overline{\widetilde{\eta}_{k}^{+}(t)\widetilde{\eta}_{k}(t)}\,\hbar\epsilon_{k}-\frac{1}{2}{\displaystyle\sum\limits_{k,}}\hbar\epsilon_{k}$
(124)

where the defining equation (46) for the gravitational modes has been used,
along with their orthogonality. Hence

$\displaystyle\left\langle\widehat{T}+\widehat{V}\right\rangle$
$\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+{\displaystyle\sum\limits_{k,}}\hbar\epsilon_{k}\,\left(\overline{\widetilde{\eta}_{k}^{+}(t)\widetilde{\eta}_{k}(t)}\,-\frac{1}{2}\right)$
(125)

The mean value for the interaction energy $\widehat{U}$ is left unchanged.

##  9 Appendix S2 - Floquet Mode Treatment

###  9.1 Position Probability Density and QCF - Floquet Modes

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

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

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 Floquet mode
functions as

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

The quantity
$\left[\overline{\widetilde{\alpha}_{k}(t)\,\widetilde{\alpha}_{l}^{+}(t)}\,-\frac{1}{2}\delta_{k,l}\right]$
is the $k,l$ element of a Hermitian matrix $H$, since
$\overline{\widetilde{\alpha}_{k}(t)\,\widetilde{\alpha}_{l}^{+}(t)}=(<\widehat{a}_{l}^{{\dagger}}$
$\widehat{a}_{k}\,>+<\widehat{a}_{k}\,\widehat{a}_{l}^{{\dagger}}>)/2$.

###  9.2 Evolution of Stochastic Phase Space Variables for Floquet Modes

By substituting the expressions in Eq (45) for
$\widetilde{\psi}(z,t),\widetilde{\psi}^{+}(z,t)$ into the Ito SFE in Eqs.
(8), (9) and using Eq (42) for the Floquet modes together with their
orthonormality property, we obtain sets of non-linear coupled equations for
the stochastic phase space variables (SPSV)

$\displaystyle\frac{\partial}{\partial t}\widetilde{\alpha}_{k}$
$\displaystyle=-i\nu_{k}\widetilde{\alpha}_{k}-i\frac{g}{\hbar}\,{\displaystyle\sum\limits_{l,m,n}}\,L_{k,l,m,n}\,\widetilde{\alpha}_{l}^{+}\,\widetilde{\alpha}_{m}\,\widetilde{\alpha}_{n}+i\frac{g}{\hbar}\,{\displaystyle\sum\limits_{l,n}}\,L_{k,l,l,n}\,\widetilde{\alpha}_{n}$
(128) $\displaystyle\frac{\partial}{\partial t}\widetilde{\alpha}_{k}^{+}$
$\displaystyle=+i\nu_{k}\widetilde{\alpha}_{k}^{+}+i\frac{g}{\hbar}\,{\displaystyle\sum\limits_{l,m,n}}\,L_{k,l,m,n}^{\ast}\,\widetilde{\alpha}_{l}\,\widetilde{\alpha}_{m}^{+}\,\widetilde{\alpha}_{n}^{+}-i\frac{g}{\hbar}\,{\displaystyle\sum\limits_{l,n}}\,L_{k,l,l,n}^{\ast}\,\widetilde{\alpha}_{n}^{+}$
(129)

where the matrix $L$ involves integrals of products of the Floquet mode
functions, defined as

$L_{k,l,m,n}={\displaystyle\int}dz\,\phi_{k}^{\ast}(z,t)\,\phi_{l}^{\ast}(z,t)\,\phi_{m}(z,t)\,\phi_{n}(z,t)$
(130)

This matrix is time dependent and periodic with period $T$.

Though non-linear, the equations for the SPSV are deterministic and can be
solved numerically if the initial values
$\widetilde{\alpha}_{k}(0),\widetilde{\alpha}_{k}^{+}(0)$ are known. These
initial values are of course stochastic and lead to the field functions given
by Eq (45) being stochastic. The distribution function for the initial values
is chosen to represent the known properties of the initial quantum state. Note
that the $\widetilde{\alpha}_{k}$ and $\widetilde{\alpha}_{k}^{+}$ SPSV do not
evolve independently. The coupled equations for the SPSV no longer involve the
periodic potential directly - this has been taken into account via the
introduction of the Floquet mode functions $\phi_{k}(z,t)$ and the Floquet
frequencies $\nu_{k}$, resulting in equations that are now focused on the
many-body quantum effects.

An alternative way of writing the non-linear coupled eqations is

$\displaystyle\frac{\partial}{\partial t}\widetilde{\alpha}_{k}$
$\displaystyle=-i\nu_{k}\widetilde{\alpha}_{k}-i\frac{g}{\hbar}{\displaystyle\sum\limits_{n}}D_{k,n}\,\widetilde{\alpha}_{n}$
(131) $\displaystyle\frac{\partial}{\partial t}\widetilde{\alpha}_{k}^{+}$
$\displaystyle=+i\nu_{k}\widetilde{\alpha}_{k}^{+}+i\frac{g}{\hbar}{\displaystyle\sum\limits_{n}}D_{k,n}^{+}\,\widetilde{\alpha}_{n}^{+}$
(132)

where

$\displaystyle D_{k,n}$
$\displaystyle={\displaystyle\int}dz\,\phi_{k}^{\ast}(z,t)\left(\widetilde{\psi}^{+}(z,t)\,\widetilde{\psi}(z,t)\,-\delta_{C}(z,z)\right)\,\phi_{n}(z,t)$
$\displaystyle D_{k,n}^{+}$
$\displaystyle={\displaystyle\int}dz\,\phi_{k}(z,t)\left(\widetilde{\psi}^{+}(z,t)\,\widetilde{\psi}(z,t)\,-\delta_{C}(z,z)\right)\,\phi_{n}^{\ast}(z,t)$
(133)

This alternative form is more convenient for numerical calculations. The non-
linearity is embodied in the matrices $D,D^{+\text{ }}.$ In this method the
stochastic fields are determined at each time point from Eq. (45), which then
can also be used to determine the quantum depletion (see Eq. (38)).

###  9.3 Initial Conditions - Floquet Modes

These will be specified via the initial stochastic field

$\widetilde{\psi}(z,0)={\displaystyle\sum\limits_{k}}\widetilde{\alpha}_{k}(0)\,\phi_{k}(z,0)\qquad\widetilde{\psi}^{+}(z,0)={\displaystyle\sum\limits_{k}}\widetilde{\alpha}_{k}^{+}(0)\,\phi_{k}^{\ast}(z,0)$
(134)

which require knowing the Floquet mode functions at $t=0$ and choosing a
stochastic distribution of the
$\widetilde{\alpha}_{k}(0),\widetilde{\alpha}_{k}^{+}(0)$ to match the quantum
state that has been prepared in the trapping potential.

##  10 Appendix S3 - Bogoliubov Theory

As foreshadowed in Sect. 2.4, in this Appendix we outline the derivation of
the Bogoliubov theory form of the Grand Canonical Hamiltonian, showing that it
is the sum of Hamiltonians for independent quantum harmonic oscilators for
each Bogoliubov mode, plus a term for the energy of the condensate mode and
some unimportant constant terms.

###  10.1 Grand Canonical Hamiltonian and Bogoliubov Approximation

The Hamiltonian $\widehat{H}$ describing the evolution of the quantum state is
given by Eq. (1), but with $V(z,t)$ replaced by $V_{trap}(z)$. If the quantum
state $\widehat{\rho}$ is invariant under the $U(1)$ symmetry group of phase
changing unitary operators $\widehat{U}(\theta)=\exp(-i\widehat{N}\,\theta)$,
it follows that $[\widehat{N},\widehat{\rho}]=0$. Hence the evolution for
$\widehat{\rho}$ can be described by the grand canonical Hamiltonian
$\widehat{K}=\widehat{H}-\mu\widehat{N}$, with $\mu$ chosen so that
$\left\langle\widehat{N}\right\rangle=N_{c}$.

In Bogoliubov theory the grand canonical Hamiltonian $\widehat{K}$ is expanded
correct to the second order in the fluctuation field and by applying the
Bogoliubov approximation in which quantum fluctuations of the condensate mode
are ignored by replacing $\widehat{c}_{0}$ with $\sqrt{N_{c}}$ giving

$\displaystyle\widehat{K}$
$\displaystyle={\displaystyle\int}dz\,\Phi_{c}(z)^{\ast}\left[-\frac{\hslash^{2}}{2m}\frac{\partial^{2}}{\partial
z^{2}}+V_{trap}(z)+\frac{1}{2}g\,n_{C}(z)-\mu\right]\,\Phi_{c}(z)$
$\displaystyle+{\displaystyle\int}dz\,\delta\widehat{\Psi}(z)\,\left[-\frac{\hslash^{2}}{2m}\frac{\partial^{2}}{\partial
z^{2}}+V_{trap}(z)+g\,n_{C}(z)-\mu\right]\,\Phi_{c}(z)^{\ast}$
$\displaystyle+{\displaystyle\int}dz\,\delta\widehat{\Psi}(z)^{{\dagger}}\,\left[-\frac{\hslash^{2}}{2m}\frac{\partial^{2}}{\partial
z^{2}}+V_{trap}(z)+g\,n_{C}(z)-\mu\right]\,\Phi_{c}(z)\,$
$\displaystyle+{\displaystyle\int}dz\,\delta\widehat{\Psi}(z)^{{\dagger}}\left[-\frac{\hslash^{2}}{2m}\frac{\partial^{2}}{\partial
z^{2}}+V_{trap}(z)+2g\,n_{C}(z)-\mu\right]\,\delta\widehat{\Psi}(z)$
$\displaystyle+{\displaystyle\int}dz\,\frac{1}{2}g\,\left[\left\\{\Phi_{c}(z)^{\ast}\delta\widehat{\Psi}(z)\right\\}^{2}+\left\\{\Phi_{c}(z)\,\delta\widehat{\Psi}(z)^{{\dagger}}\right\\}^{2}\right]$
(135)

This is equivalent to writing the field operator as
$\hat{\Psi}(z)=\Phi_{c}(z)+\delta\widehat{\Psi}(z)$, so the condensate field
term is replaced by a non-operator field
$\sqrt{N_{c}}\psi_{c}(z)=\Phi_{c}(z)$. Note however that we still will require
the condensate mode to be orthogonal to the modes associated with the
fluctuation field, and the commutation rules in Eq. (67) to still apply. In
the expression for $\widehat{K}$ the first line gives the zero order
contribution as a constant term, the next two lines the first order
contribution and the last two lines the second order contribution. The boson
number density associated with the condensate is given by
$n_{C}(z)=\Phi_{c}(z)^{\ast}\Phi_{c}(z)$. Note that a further approximation
has been made - there are terms involving
$\,(\delta\widehat{\Psi}(z)^{{\dagger}})^{2}\,\delta\widehat{\Psi}(z)^{2},\,(\delta\widehat{\Psi}(z)^{{\dagger}})^{2}\,\delta\widehat{\Psi}(z),\,\delta\widehat{\Psi}(z)^{{\dagger}}\,\delta\widehat{\Psi}(z)^{2}$
that have been discarded.

Since $\Phi_{c}(z)$ satisfies Eq. (60), the linear terms in the grand
canonical Hamiltonian are zero. The constant term has no dynamical effect. The
second order term that remains may be diagonalised via the Bogoliubov
transformation, which is an example of a linear canonical transformation (see
Ref. [22], Sections 6.1, 6.4.3) in which commutation rules for the mode
operators are preserved.

###  10.2 Bogoliubov Hamiltonian

On substituting for the fluctuation field in terms of the Bogoliubov mode
operators using Eq. (68) and with $u_{k}(z),v_{k}(z)$ satisfying the
generalised BDG equations (73) , the grand canonical Hamiltonian given by Eq.
(135) can be expressed as the sum of Hamiltonians for independent quantum
harmonic oscilators for each Bogoliubov mode, plus a term for the energy of
the condensate mode and some unimportant constant terms [32]. Thus

$\displaystyle\widehat{K}$
$\displaystyle={\displaystyle\int}dz\,\Phi_{c}(z)^{\ast}\left[-\frac{\hslash^{2}}{2m}\frac{\partial^{2}}{\partial
z^{2}}+V_{trap}(z)+\frac{1}{2}g\,n_{C}(z)-\mu\right]\,\Phi_{c}(z)$
$\displaystyle+{\displaystyle\sum\limits_{k\neq
0}}\left[-{\displaystyle\int}dz\,|\,V_{k}(z)|^{2}\right]$
$\displaystyle+{\displaystyle\sum\limits_{k\neq
0}}\hslash\omega_{k}\,\widehat{b}_{k}^{{\dagger}}\,\widehat{b}_{k}$ (136)

In deriving Eq (136) the Hermitiancy properties of $\mathcal{L}$, the reality
of the $\omega_{k}$ , the commutation rules for the
$\widehat{b}_{k},\widehat{b}_{k}^{{\dagger}}$, the biorthogonality conditions
(72) and the orthogonality conditions (69) are all used, in addition to the
generalised BDG equations (73).

##  11 Appendix S4 - Alternative Initial States for Condensate

One simple idea would be to assume the quantum state is a pure state given by
a single Glauber coherent state $\left|\gamma 0\right\rangle_{c}$ for the
condensate mode with amplitude $\gamma_{0}$, and with each Bogoliubov mode $k$
in its vacuum state $\left|0\right\rangle_{k}$. Thus

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

with

$\displaystyle\left|\Phi_{\gamma 0}\right\rangle$ $\displaystyle=\left|\gamma
0\right\rangle_{c}\left|0\right\rangle_{1}\left|0\right\rangle_{2}....\left|0\right\rangle_{k}....$
$\displaystyle=\left|\gamma
0\right\rangle_{c}\times{\textstyle\prod\limits_{k\neq
0}}\,\left|0\right\rangle_{k}$ (138)

For such a state we would have, if $\gamma 0\,=\sqrt{N_{c}}$

$\displaystyle\hat{\Psi}(z)\,\left|\Phi_{\gamma 0}\right\rangle$
$\displaystyle=\sqrt{N_{c}}\,\psi_{c}(z)\,\left|\Phi_{\gamma
0}\right\rangle=\Phi_{c}(z)\,\left|\Phi_{\gamma 0}\right\rangle$
$\displaystyle\left\langle\hat{\Psi}(z)\right\rangle$
$\displaystyle=Tr\,\hat{\Psi}(z)\,\widehat{\rho}(0)=\Phi_{c}(z)\,$
$\displaystyle\left\langle\hat{\Psi}(z)^{{\dagger}}\,\hat{\Psi}(z)\right\rangle$
$\displaystyle=Tr\,\hat{\Psi}(z)^{{\dagger}}\,\hat{\Psi}(z)\,\widehat{\rho}(0)=\Phi_{c}(z)^{\ast}\Phi_{c}(z)=n_{C}(z)$
$\displaystyle\left\langle\widehat{N}\right\rangle$ $\displaystyle=\int
dz\,\left\langle\hat{\Psi}(z)^{{\dagger}}\,\hat{\Psi}(z)\right\rangle=N_{c}$
(139)

using $\widehat{c}_{0}\left|\gamma 0\right\rangle_{c}=\gamma 0\left|\gamma
0\right\rangle_{c}$ for Glauber coherent states. This state does obviously
resemble a BEC with all $N_{c}$ bosons having the same wave function
$\Phi_{c}(z)$. However, it is not invariant under phase change
transformations, so this is inconsistent with the requirement of
$\widehat{\rho}$ being phase invariant in order to treat evolution via the
grand canonical Hamiltonian $\widehat{K}$.

Nevertheless, a simple quantum state that is phase invariant can easily be
constructed as a mixed state based on the $\left|\Phi_{\gamma
0}\right\rangle$. Writing $\gamma 0\,=\sqrt{N_{c}}\exp(i\phi_{c})$ we now
consider the mixed state given by

$\widehat{\rho}(0)=\int_{0}^{2\pi}\frac{d\phi_{c}}{2\pi}\left|\Phi_{\gamma
0}\right\rangle\left\langle\Phi_{\gamma 0}\right|$ (140)

We then find that

$\displaystyle\hat{\Psi}(z)\,\left|\Phi_{\gamma 0}\right\rangle$
$\displaystyle=\sqrt{N_{c}}\exp(i\phi_{c})\,\psi_{c}(z)\,\left|\Phi_{\gamma
0}\right\rangle=\exp(i\phi_{c})\,\Phi_{c}(z)\,\left|\Phi_{\gamma
0}\right\rangle$ $\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)=\Phi_{c}(z)^{\ast}\Phi_{c}(z)$ (141)

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. This mixed state is
therefore one good description of the initial BEC.

##  12 Appendix S5 - General Initial States

For the state given by Eq (140) the QCF for normally ordered products of the
field operators are given by

$\displaystyle\left\langle\hat{\Psi}(z)\right\rangle$
$\displaystyle=\left\langle\hat{\Psi}(z^{\\#})^{{\dagger}}\right\rangle=0$
(142)
$\displaystyle\left\langle\hat{\Psi}(z^{\\#})^{{\dagger}}\,\hat{\Psi}(z)\right\rangle$
$\displaystyle=\Phi_{c}(z^{\\#})^{\ast}\Phi_{c}(z)$
$\displaystyle\left\langle\hat{\Psi}(z_{1})\,\hat{\Psi}(z_{2})\right\rangle$
$\displaystyle=\left\langle\hat{\Psi}(z_{1}^{\\#})^{{\dagger}}\,\hat{\Psi}(z_{2}^{\\#})^{{\dagger}}\right\rangle=0$
(143) $\displaystyle.....$
$\displaystyle\left\langle\hat{\Psi}(z_{n}^{\\#})^{{\dagger}}\,.....\hat{\Psi}(z_{1}^{\\#})^{{\dagger}}\,\hat{\Psi}(z_{1}).....\hat{\Psi}(z_{n})\right\rangle$
$\displaystyle=\Phi_{c}(z_{n}^{\\#})^{\ast}.....\Phi_{c}(z_{1}^{\\#})^{\ast}\Phi_{c}(z_{1}).....\Phi_{c}(z_{n})$
(144)
$\displaystyle\left\langle\hat{\Psi}(z_{n}^{\\#})^{{\dagger}}\,.....\hat{\Psi}(z_{1}^{\\#})^{{\dagger}}\,\hat{\Psi}(z_{1})......\hat{\Psi}(z_{m})\right\rangle$
$\displaystyle=0\qquad n\neq m$ (145)

Note that the normally ordered QCF is zero unless there are the same numbers
of $\hat{\Psi}(z)$ and $\hat{\Psi}(z)^{{\dagger}}$. When the numbers are the
same, the QCF is determined from the condensate wave function.

Normally ordered QCF are related to symmetrically ordered QCF via expressions
such as $\\{\hat{\Psi}^{{\dagger}}(\ z^{\\#})\hat{\Psi}(\
z)\\}=\hat{\Psi}^{{\dagger}}(\ z^{\\#})\hat{\Psi}(\
z)+\frac{1}{2}\delta_{C}(z,z^{\\#})$. This is a particular case of Wick’s
theorem which states that a product of field operators such as
$\hat{\Psi}(z_{n}^{\\#})^{{\dagger}}\,.....\hat{\Psi}(z_{1}^{\\#})^{{\dagger}}\,\hat{\Psi}(z_{1})......\hat{\Psi}(z_{m})$
may be written as

$\displaystyle\widehat{A}_{1}\widehat{A}_{2}.....\widehat{A}_{n}$
$\displaystyle=N[\widehat{A}_{1}\widehat{A}_{2}....\widehat{A}_{n}]+{\textstyle\sum\limits_{(i,j)}}N[\widehat{A}_{1}\widehat{A}_{2}..\overbrace{\widehat{A}_{i}...\widehat{A}_{j}}..\widehat{A}_{n}]+{\textstyle\sum\limits_{(i,j)(k,l)}}N[\widehat{A}_{1}\widehat{A}_{2}..\overbrace{\widehat{A}_{i}...\widehat{A}_{j}}...\overbrace{\widehat{A}_{k}...\widehat{A}_{l}}....\widehat{A}_{n}]+...$
(146)

where $N[...]$ is the normally ordered form of the quantity in square
brackets, and $\overbrace{\widehat{A}_{i}\widehat{A}_{j}}$ is the so-called
contraction of the pair of operators $\widehat{A}_{i}\widehat{A}_{j}$ which is
defined by
$\overbrace{\widehat{A}_{i}\widehat{A}_{j}}=\widehat{A}_{i}\widehat{A}_{j}-N[\widehat{A}_{i}\widehat{A}_{j}]$.
In the second and subsequent terms on the right side, the contractions (which
are c-numbers ) are removed as multiplying factors to the remaining normalised
ordered expression. In our case the relevant contractions are

$\displaystyle\overbrace{\hat{\Psi}(z_{1}^{\\#})^{{\dagger}}\hat{\Psi}(z_{2}^{\\#})^{{\dagger}}}$
$\displaystyle=0\qquad\qquad\overbrace{\hat{\Psi}(z_{1})\hat{\Psi}(z_{2})}=0$
(147)
$\displaystyle\overbrace{\hat{\Psi}(z_{1}^{\\#})^{{\dagger}}\hat{\Psi}(z_{1})}$
$\displaystyle=0\qquad\qquad\overbrace{\hat{\Psi}(z_{1})\hat{\Psi}(z_{1}^{\\#})^{{\dagger}}}=\delta(z_{1}-z_{1}^{\\#})$

Hence we have

$\displaystyle\overline{\widetilde{\psi}(z)}$
$\displaystyle=\left\langle\left\\{\hat{\Psi}(z)\right\\}\right\rangle=\left\langle\hat{\Psi}(z)\right\rangle=0$
$\displaystyle\overline{\widetilde{\psi}^{+}(z)}$
$\displaystyle=\left\langle\left\\{\hat{\Psi}(z)^{{\dagger}}\right\\}\right\rangle=\left\langle\hat{\Psi}(z)^{{\dagger}}\right\rangle=0$
(148)

for first order QCF, and

$\displaystyle\overline{\widetilde{\psi}(z_{1})\widetilde{\psi}(z_{2})}$
$\displaystyle=\left\langle\left\\{\hat{\Psi}(z_{1})\hat{\Psi}(z_{2})\right\\}\right\rangle=\left\langle\hat{\Psi}(z_{1})\hat{\Psi}(z_{2})\right\rangle=0$
$\displaystyle\overline{\widetilde{\psi}^{+}(z_{1}^{\\#})\widetilde{\psi}^{+}(z_{2}^{\\#})}$
$\displaystyle=\left\langle\left\\{\hat{\Psi}(z_{1}^{\\#})^{{\dagger}}\hat{\Psi}(z_{2}^{\\#})^{{\dagger}}\right\\}\right\rangle=\left\langle\hat{\Psi}(z_{1}^{\\#})^{{\dagger}}\hat{\Psi}(z_{2}^{\\#})^{{\dagger}}\right\rangle=0$
$\displaystyle\overline{\widetilde{\psi}^{+}(z_{1}^{\\#})\widetilde{\psi}^{+}(z_{1})}$
$\displaystyle=\left\langle\left\\{\hat{\Psi}(z_{1}^{\\#})^{{\dagger}}\hat{\Psi}(z_{1})\right\\}\right\rangle=\left\langle\hat{\Psi}(z_{1}^{\\#})^{{\dagger}}\hat{\Psi}(z_{1})\right\rangle+\delta(z_{1}-z_{1}^{\\#})=\Phi_{c}(z_{1}^{\\#})^{\ast}\Phi_{c}(z_{1})+\delta(z_{1}-z_{1}^{\\#})$
(149)

for second order QCF.

By substituting for the stochastic fields from Eq (82) (reverting to the
$u_{k}$ and $v_{k}$ notation for the Bogoliubov modes) and making use of the
orthogonality properties (69) of the condensate and Bogoliubov modes, we can
show from the first order QCF results (148) that

$\displaystyle\overline{\widetilde{\gamma}_{0}(0)}$
$\displaystyle=\overline{\widetilde{\gamma}_{0}^{+}(0)}=0$
$\displaystyle\overline{\,\widetilde{\beta}_{k}(0)}$
$\displaystyle=\overline{\,\widetilde{\beta}_{k}^{+}(0)}=0$ (150)

showing that the stochastic averages of the initial stochastic amplitudes for
the condensate mode and the Bogoliubov modes are all zero.

Similarly, by making use of both the orthogonality properties (69) of the
condensate and Bogoliubov modes along with the biorthogonality conditions (72)
for the Bogoliubov modes, we can show from the second order QCF results (149)
that

$\displaystyle\overline{\widetilde{\gamma}_{0}^{+}(0)\widetilde{\gamma}_{0}(0)}$
$\displaystyle=N_{c}+\frac{1}{2}\qquad\overline{\widetilde{\gamma}_{0}(0)\widetilde{\gamma}_{0}(0)}=0\qquad\overline{\widetilde{\gamma}_{0}^{+}(0)\widetilde{\gamma}_{0}^{+}(0)}=0$
(151)
$\displaystyle\overline{\widetilde{\gamma}_{0}(0)\widetilde{\beta}_{k}(0)}$
$\displaystyle=\overline{\widetilde{\gamma}_{0}(0)\widetilde{\beta}_{k}^{+}(0)}=\overline{\widetilde{\gamma}_{0}^{+}(0)\widetilde{\beta}_{k}(0)}=\overline{\widetilde{\gamma}_{0}^{+}(0)\widetilde{\beta}_{k}^{+}(0)}=0$
(152)

This shows that the quantity
$\overline{\widetilde{\gamma}_{0}^{+}(0)\widetilde{\gamma}_{0}(0)}$ is never
less than one half - reflecting the quantum nature of the condensate mode, and
increases linearly with the number of bosons in this mode. The second set of
results reflects the lack of initial correlation between the condensate and
Bogoliubov modes.

The second order QCF results (149) also could be used to find expressions for
stochastic averages involving pairs of Bogoliubov modes such as
$\,\overline{\widetilde{\beta}_{k}^{+}(0)\widetilde{\beta}_{m}(0)}$. However,
this is more easily accomplished by considering the modes separately.

##  13 Appendix S6 - Wigner Distribution Function

###  13.1 Factorisation of Wigner Distribution Function

We first show that the Wigner distribution function (see Sections 7.2, 7.3 in
Ref [22] for basic definiions of Wigner distribution functions) for condensate
and non-condensate modes factorises into a Wigner distribution function
$W_{C}(\gamma_{0},\gamma_{0}^{+})$ for the condensate mode and a Wigner
distribution function $W_{NC}(\mathbf{\gamma},\mathbf{\gamma}^{+})$ for the
non-condensate modes, where $\gamma_{0},\gamma_{0}^{+}$ are the phase space
variables for the condensate mode representing
$\widehat{c}_{0},\widehat{c}_{0}^{{\dagger}}$ and where
$\mathbf{\gamma\equiv\\{\gamma}_{1},...,\gamma_{i},...\\}$ and
$\mathbf{\gamma}^{+}\mathbf{\equiv\\{\gamma}_{1}^{+},...,\gamma_{i}^{+},...\\}$
are the phase space variables for the standard non-condensate modes
$\widehat{c}_{i},\widehat{c}_{i}^{{\dagger}}$ $(i\neq 0)$. We begin with the
characteristic function $\chi_{NC}(\xi_{0,}\xi_{0}^{+},\mathbf{\xi,\xi}^{+})$
for the condensate mode $\widehat{c}_{0},\widehat{c}_{0}^{{\dagger}}$ and
standard non-condensate modes $\widehat{c}_{i},\widehat{c}_{i}^{{\dagger}}$
$(i\neq 0)$, which is defined by

$\chi_{W}(\xi_{0,}\xi_{0}^{+},\mathbf{\xi,\xi}^{+})=Tr\,\exp((i\\{\widehat{c}_{0}\times\xi_{0}^{+}+\xi_{0}\times\widehat{c}_{0}^{{\dagger}}+\widehat{\mathbf{c}}\cdot\mathbf{\xi}^{+}+\mathbf{\xi\cdot}\widehat{\mathbf{c}}^{{\dagger}}\\})\,\widehat{\rho})$
(153)

where
$\widehat{\mathbf{c}}\equiv\\{\widehat{c}_{1},...,\widehat{c}_{i},....\\}$,
$\widehat{\mathbf{c}}^{{\dagger}}\equiv\\{\widehat{c}_{1}^{{\dagger}},...,\widehat{c}_{i}^{{\dagger}},....\\}$,
$\mathbf{\xi\equiv\\{\xi}_{1},...,\xi_{i},...\\}$ and
$\mathbf{\xi}^{+}\mathbf{\equiv\\{\xi}_{1}^{+},...,\xi_{i}^{+},...\\}$. This
is related to the Wigner distribution function by

$\chi_{W}(\xi_{0,}\xi_{0}^{+},\mathbf{\xi,\xi}^{+})={\textstyle\int}d^{2}\gamma_{0}\,d^{2}\gamma_{0}^{+}\,d^{2}\mathbf{\gamma\,}d^{2}\mathbf{\gamma}^{+}\mathbf{\,}\exp(i\\{\gamma_{0}\,\xi_{0}^{+}+\xi_{0}\,\gamma_{0}^{+}+\mathbf{\gamma}\cdot\mathbf{\xi}^{+}+\mathbf{\xi\cdot\gamma}^{+}\\})\,W(\gamma_{0},\gamma_{0}^{+},\mathbf{\gamma},\mathbf{\gamma}^{+})$
(154)

Now since $\widehat{\rho}=\widehat{\rho}_{C}\otimes\widehat{\rho}_{NC}$ we see
that the characteristic function factorises as

$\displaystyle\chi_{W}(\xi_{0,}\xi_{0}^{+},\mathbf{\xi,\xi}^{+})$
$\displaystyle=Tr_{C}\,\exp((i\\{\widehat{c}_{0}\times\xi_{0}^{+}+\xi_{0}\times\widehat{c}_{0}^{{\dagger}}\\})\,\widehat{\rho}_{C})\times
Tr_{NC}\,\exp((i\\{\widehat{\mathbf{c}}\cdot\mathbf{\xi}^{+}+\mathbf{\xi\cdot}\widehat{\mathbf{c}}^{{\dagger}}\\})\,\widehat{\rho}_{NC})$
$\displaystyle=\chi_{C}(\xi_{0,}\xi_{0}^{+})\times\chi_{NC}(\mathbf{\xi,\xi}^{+})$
(155)

Hence we can write

$\displaystyle\chi_{C}(\xi_{0,}\xi_{0}^{+})$
$\displaystyle={\textstyle\int}d^{2}\gamma_{0}\,d^{2}\gamma_{0}^{+}\,\exp(i\\{\gamma_{0}\,\xi_{0}^{+}+\xi_{0}\,\gamma_{0}^{+}\\})\,W_{C}(\gamma_{0},\gamma_{0}^{+})$
(156) $\displaystyle\chi_{NC}(\mathbf{\xi,\xi}^{+})$
$\displaystyle={\textstyle\int}d^{2}\mathbf{\gamma\,}d^{2}\mathbf{\gamma}^{+}\mathbf{\,}\exp(i\\{\mathbf{\gamma}\cdot\mathbf{\xi}^{+}+\mathbf{\xi\cdot\gamma}^{+}\\})\,W_{NC}(\mathbf{\gamma},\mathbf{\gamma}^{+}$
(157)

This then defines the Wigner distribution functions for the condensate mode
and for the non-condensate modes.

###  13.2 Relation between Standard and Bogoliubov Non-Condensate Modes

By equating expressions (63) and (68) for the fluctuation field operator in
terms of standard modes $\psi_{i}(z)$ and Bogoliubov modes $u_{k}(z),v_{k}(z)$
and using the orthogonality of the standard modes we follow the approach of
Morgan [25] and express the standard non-condensate mode annihilation,
creation operators $\widehat{c}_{i},\widehat{c}_{i}^{{\dagger}}$ in terms of
the Bogoliubov operators $\widehat{b}_{k},\widehat{b}_{k}^{{\dagger}}$. In
column matrix form this is

$\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]$ (158)

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}$ (159)

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.

Using the bosonic commutation rules for
$\widehat{c}_{i},\widehat{c}_{i}^{{\dagger}}$ and for
$\widehat{b}_{k},\widehat{b}_{k}^{{\dagger}}$ leads to the following matrix
equations

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

It then follows from the last equations that

$\left[\begin{array}[c]{cc}U&-V\\\
-V^{\ast}&U^{\ast}\end{array}\right]\times\left[\begin{array}[c]{cc}U^{{\dagger}}&V^{T}\\\
V^{{\dagger}}&U^{T}\end{array}\right]=\left[\begin{array}[c]{cc}E&0\\\
0&E\end{array}\right]$ (161)

We can then use the last equation to express the Bogoliubov operators
$\widehat{b}_{k},\widehat{b}_{k}^{{\dagger}}$ in terms of the standard mode
annihilation, creation operators
$\widehat{c}_{i},\widehat{c}_{i}^{{\dagger}}$. This gives

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

and finally we can also derive a further result involving the $U,V$ matrices

$\left[\begin{array}[c]{cc}U&-V\\\
-V^{\ast}&U^{\ast}\end{array}\right]^{({\dagger})}\times\left[\begin{array}[c]{cc}E&0\\\
0&-E\end{array}\right]\times\left[\begin{array}[c]{cc}U&-V\\\
-V^{\ast}&U^{\ast}\end{array}\right]=\left[\begin{array}[c]{cc}E&0\\\
0&-E\end{array}\right]$ (163)

The last result involves evaluating the matrix elements for
$(U^{{\dagger}}U-V^{T}V^{\ast})$ and $(-V^{T}U+U^{T}V^{\ast})$, using Eqs.
(159) for the matrix elements $U_{i,k}$ and $V_{i,k}$, together with the
completeness relation
${\displaystyle\sum\limits_{i}}\psi_{i}(z)\,\psi_{i}^{\ast}(z^{\\#})=\delta(z-z^{\\#})$
for standard mode functions as well as the biorthogonality conditions (72) and
the condensate mode orthogonality conditions (69).

###  13.3 Wigner Distribution Function for Bogoliubov Modes

To show that a Wigner distribution for the non-condensate modes exists in
terms of phase variables $\mathbf{\beta}\equiv\\{\beta_{1},...\beta_{k},..)$
and $\mathbf{\beta}^{+}\equiv\\{\beta_{1}^{+},...\beta_{k}^{+},..)$ for the
Bogoliubov modes, we begin with the characteristic function
$\chi_{NC}(\mathbf{\xi,\xi}^{+})$ for the standard modes
$\widehat{c}_{i},\widehat{c}_{i}^{{\dagger}}$ introduced previously

$\chi_{NC}(\mathbf{\xi,\xi}^{+})=Tr_{NC}(\exp(i\\{\widehat{\mathbf{c}}\cdot\mathbf{\xi}^{+}+\mathbf{\xi\cdot}\widehat{\mathbf{c}}^{{\dagger}}\\})\,\widehat{\rho}_{NC})$
(164)

where $\mathbf{\xi\equiv\\{\xi}_{1},...,\xi_{i},...\\}$ and
$\mathbf{\xi}^{+}\mathbf{\equiv\\{\xi}_{1}^{+},...,\xi_{i}^{+},...\\}$. This
is related to the Wigner distribution function by

$\chi_{NC}(\mathbf{\xi,\xi}^{+})={\textstyle\int}d^{2}\mathbf{\gamma\,}d^{2}\mathbf{\gamma}^{+}\mathbf{\,}\exp(i\\{\mathbf{\gamma}\cdot\mathbf{\xi}^{+}+\mathbf{\xi\cdot\gamma}^{+}\\})\,W_{NC}(\mathbf{\gamma},\mathbf{\gamma}^{+})$
(165)

where $\mathbf{\gamma\equiv\\{\gamma}_{1},...,\gamma_{i},...\\}$ and
$\mathbf{\gamma}^{+}\mathbf{\equiv\\{\gamma}_{1}^{+},...,\gamma_{i}^{+},...\\}$
are the phase space variables for the standard modes.

Now the phase space variables $\mathbf{\gamma},\mathbf{\gamma}^{+}$ associated
with the standard modes will be related to the phase space variables
$\mathbf{\beta},\mathbf{\beta}^{+}$ associated with the Bogoliubov modes via
the same form (158) that applies for the mode annihilation, creation
operators. Thus

$\left[\begin{array}[c]{c}\mathbf{\gamma}\\\
\mathbf{\gamma}^{+}\end{array}\right]=\left[\begin{array}[c]{cc}U&-V\\\
-V^{\ast}&U^{\ast}\end{array}\right]\times\left[\begin{array}[c]{c}\mathbf{\beta}\\\
\mathbf{\beta}^{+}\end{array}\right]$ (166)

The same applies to the stochastic phase space variables. Hence we have

$\displaystyle\left[\begin{array}[c]{c}\widetilde{\mathbf{\gamma}}\\\
\widetilde{\mathbf{\gamma}}^{+}\end{array}\right]$
$\displaystyle=\left[\begin{array}[c]{cc}U&-V\\\
-V^{\ast}&U^{\ast}\end{array}\right]\times\left[\begin{array}[c]{c}\widetilde{\mathbf{\beta}}\\\
\widetilde{\mathbf{\beta}}^{+}\end{array}\right]$ (173)
$\displaystyle\left[\begin{array}[c]{c}\widetilde{\mathbf{\beta}}\\\
\widetilde{\mathbf{\beta}}^{+}\end{array}\right]$
$\displaystyle=\left[\begin{array}[c]{cc}U^{{\dagger}}&V^{T}\\\
V^{{\dagger}}&U^{T}\end{array}\right]\times\left[\begin{array}[c]{c}\widetilde{\mathbf{\gamma}}\\\
\widetilde{\mathbf{\gamma}}^{+}\end{array}\right]$ (180)

where we have also written down the inverse relation, based on (162).

Using the result (158) we can write
$\exp(i\\{\mathbf{\gamma}\cdot\mathbf{\xi}^{+}+\mathbf{\xi\cdot\gamma}^{+}\\})=\exp(i\\{\mathbf{\beta}\cdot\mathbf{\eta}^{+}+\mathbf{\eta\cdot\beta}^{+}\\})$,
where

$\left[\begin{array}[c]{c}\mathbf{\eta}\\\
\mathbf{\eta}^{+}\end{array}\right]=\left[\begin{array}[c]{cc}U^{{\dagger}}&-V^{T}\\\
-V^{{\dagger}}&U^{T}\end{array}\right]\times\left[\begin{array}[c]{c}\mathbf{\xi}\\\
\mathbf{\xi}^{+}\end{array}\right]$ (181)

As the differentials transform as

$d^{2}\mathbf{\gamma\,}d^{2}\mathbf{\gamma}^{+}\mathbf{\,=}\left|\left[\begin{array}[c]{cc}U&-V\\\
-V^{\ast}&U^{\ast}\end{array}\right]\right|\times
d^{2}\beta\mathbf{\,}d^{2}\mathbf{\beta}^{+}\mathbf{\,}$ (182)

we can then use the result (163) to show that the determinant in (182) is
merely a constant of magnitude one.

Since the standard phase space and characteristic function variables are
linear functions of the corresponding Bogoliubov variables and also have
$(\exp(i\\{\widehat{\mathbf{b}}\cdot\mathbf{\eta}^{+}+\eta\mathbf{\cdot}\widehat{\mathbf{b}}^{{\dagger}}\\})=(\exp(i\\{\widehat{\mathbf{c}}\cdot\mathbf{\xi}^{+}+\mathbf{\xi\cdot}\widehat{\mathbf{c}}^{{\dagger}}\\})$
we thus can replace Eq. (165) by

$\displaystyle\chi_{NC}(\mathbf{\eta,\eta}^{+})$
$\displaystyle=Tr_{NC}(\exp(i\\{\widehat{\mathbf{b}}\cdot\mathbf{\eta}^{+}+\eta\mathbf{\cdot}\widehat{\mathbf{b}}^{{\dagger}}\\})\,\widehat{\rho}_{NC})$
(183)
$\displaystyle={\textstyle\int}d^{2}\beta\mathbf{\,}d^{2}\mathbf{\beta}^{+}\mathbf{\,}\exp(i\\{\mathbf{\beta}\cdot\mathbf{\eta}^{+}+\mathbf{\eta\cdot\beta}^{+}\\}\,W_{NC}(\mathbf{\beta},\mathbf{\beta}^{+})$
(184)

showing that a Wigner distribution function can be defined in terms of
Bogoliubov modes, even though these do not satisfy standard orthogonality
conditions.

###  13.4 Non-Condensate Modes and Squeezed Vacuum State

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 that the vacuum state for Bogoliubov modes is equivalent to a
squeezed vacuum state for the standard non-condensate modes.

For the first order QCF we have from Eq. (173)

$\overline{\left[\begin{array}[c]{c}\widetilde{\mathbf{\gamma}}\\\
\widetilde{\mathbf{\gamma}}^{+}\end{array}\right]}=\left[\begin{array}[c]{cc}U&-V\\\
-V^{\ast}&U^{\ast}\end{array}\right]\times\overline{\left[\begin{array}[c]{c}\widetilde{\mathbf{\beta}}\\\
\widetilde{\mathbf{\beta}}^{+}\end{array}\right]}$

so that as
$\overline{\widetilde{\beta}_{k}}=\overline{\widetilde{\beta}_{k}^{+}}=0$ it
follows that
$\overline{\widetilde{\gamma}_{i}}=\overline{\widetilde{\gamma}_{i}^{+}}=0$.
There is no difference between the standard and Bogoliubov non-condensate
modes in this regard. However, this is not the case for second order QCF. For
the second order QCF we see that from Eq. (173) that

$\displaystyle\overline{\left[\begin{array}[c]{c}\widetilde{\mathbf{\gamma}}^{+}\\\
\widetilde{\mathbf{\gamma}}\end{array}\right]\times\left[\begin{array}[c]{cc}\widetilde{\mathbf{\gamma}}&\widetilde{\mathbf{\gamma}}^{+}\end{array}\right]}$
(188) $\displaystyle=\left[\begin{array}[c]{cc}U^{\ast}&-V^{\ast}\\\
-V&U\end{array}\right]\times\overline{\left[\begin{array}[c]{c}\widetilde{\mathbf{\beta}}^{+}\\\
\widetilde{\mathbf{\beta}}\end{array}\right]\times\left[\begin{array}[c]{cc}\widetilde{\mathbf{\beta}}&\widetilde{\mathbf{\beta}}^{+}\end{array}\right]}\times\left[\begin{array}[c]{cc}U^{T}&-V^{{\dagger}}\\\
-V^{T}&U^{{\dagger}}\end{array}\right]$ (196)

where the stochastic averages of both sides have been taken. But from (103) we
have

$\overline{\left[\begin{array}[c]{c}\widetilde{\mathbf{\beta}}^{+}\\\
\widetilde{\mathbf{\beta}}\end{array}\right]\times\left[\begin{array}[c]{cc}\widetilde{\mathbf{\beta}}&\widetilde{\mathbf{\beta}}^{+}\end{array}\right]}=\frac{1}{2}\left[\begin{array}[c]{cc}E&0\\\
0&E\end{array}\right]$ (197)

so that

$\overline{\left[\begin{array}[c]{c}\widetilde{\mathbf{\gamma}}^{+}\\\
\widetilde{\mathbf{\gamma}}\end{array}\right]\times\left[\begin{array}[c]{cc}\widetilde{\mathbf{\gamma}}&\widetilde{\mathbf{\gamma}}^{+}\end{array}\right]}=\frac{1}{2}\left[\begin{array}[c]{cc}(UU^{{\dagger}}+VV^{{\dagger}})^{\ast}&(-VU^{T}-UV^{T})^{\ast}\\\
-VU^{T}-UV^{T}&(UU^{{\dagger}}+VV^{{\dagger}})\end{array}\right]$ (198)

Evaluating these sub-matrices using (159) and (70) gives

$\displaystyle\frac{1}{2}(UU^{{\dagger}}+VV^{{\dagger}})_{i,j}$
$\displaystyle=\frac{1}{2}\delta_{i,j}+(VV^{{\dagger}})_{i,j}$
$\displaystyle-\frac{1}{2}(VU^{T}+UV^{T})_{i,j}$
$\displaystyle=-(VU^{T})_{i,j}$ (199)

so that
$\overline{\widetilde{\gamma}_{i}^{+}\widetilde{\gamma}_{i}}=\frac{1}{2}+(VV^{{\dagger}})_{i,i}$,
which is always greater than $\frac{1}{2}$. This shows that the state with all
Bogoliubov modes in the vacuum state is not a state where the standard non-
condensate modes are in the vacuum state. Also
$\overline{\widetilde{\gamma}_{i}\widetilde{\gamma}_{i}}=-(VU^{T})_{i,i}$
which is no longer zero as for the Bogoliubov modes. Described in terms of the
standard non-condensate modes, the non-condensate state is actually a multi-
mode squeezed vacuum.

##  14 Appendix S7 - No Driving Case

To confirm that both driving and interactions must both be present for DTC
creation, we consider Figs. 14, 15 for the PPD and OBP for the case of
interactions $(gN=-0.01)$ but no driving, and Figs. 4, 9 for the case of no
interaction $(gN=0)$ but with driving. In the latter case the PPD shows a
mixture of $T,2T$ periodicity corresponding to transfer of bosons back and
forth between the two Wannier modes, so no simple $2T$ periodicity occurs. In
the former case, the PPD shows an irregular behaviour after a transient
interval where a $2T$ periodicity (associated with t${}_{bounce}=2T$) is
initially seen, but which rapidly disappears. The corresponding OBP and its FT
shows that there is no regular periodicity. Unlike cases where both driving
and interactions are present and a DTC occurs, a DTC is not present in
situations where one factor is absent. Of course just having both factors
present does not guarantee DTC behaviour. The interaction needs to be
sufficiently strong to allow DTC behaviour.

Figure 14: PPD for harmonic trap initial conditions and static mirror case as
a function of $z$ and $t$ for $gN=-0.1$ and $N=600$. The blue solid (red
dashed) curves are calculated using the TWA (GPE) approach. Figure 15: OBP and
the corresponding FT for harmonic trap initial conditions and static mirror
case with $gN=-0.1$ and $N=600$. The blue circle (red cross) symbols show the
TWA (GPE) results.