Many-Body Effects and Quantum Fluctuations for Discrete Time Crystals in Bose-Einstein Condensates

We present a fully comprehensive multi-mode quantum treatment based on the truncated Wigner approximation (TWA) to study many-body effects and quantum fluctuations on the formation of a discrete time crystal (DTC) in a Bose-Einstein condensate (BEC) bouncing resonantly on an atom mirror, driven at period T. 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. 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 T, 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 600 atoms at times corresponding to 2000 T, 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 of the DTC does not increase significantly for times out to at least 2000 mirror oscillations, so TWA theory predicts that thermalisation is absent. Finally, we find that the dynamical behaviour is similar for attractive or repulsive boson-boson interactions, and that a stable DTC based on repulsive interactions can be created.

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 timedependent 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.

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 meanfield 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 manybody 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 phasespace 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 ≤ |gN | ≤ 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.

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].

Hamiltonian
The Hamiltonian for the 1D model for the periodically driven BEC is given in terms of the field operatorsΨ(z),Ψ(z) † for the annihilation, creation of a bosonic atom of mass m at position z as [17] and is the sum of kinetic energy, time-dependent potential energy and interaction energy terms. The potential energy contains a driving term at frequency ω = 2π/T and is given by 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 λ. 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 ω ⊥ a s -which is proportional to the s-wave scattering length a s and the oscillation frequency ω ⊥ for the BEC atoms in a transverse trap -may be tuned via a Feshbach resonance.

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.

Functional Fokker-Planck Equation -W + Distribution
The functional Fokker-Planck equation (FFPE) for the W + Wigner distribution functional W [ψ(z)] which represents the quantum density operator ρ 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 K, V , U . A derivation for 3D is presented in Ref [22] (see Section 15.1.6 and Appendix I). The field operatorsΨ(z),Ψ(z) † are represented by two unrelated time-independent c-number fields ψ(z) and ψ + (z), respectively. For short, ψ(z) ≡ ψ(z), ψ + (z), and W [ψ(z)] is time dependent.
The kinetic energy term is and only contributes to the drift term in the FFPE. The potential energy term is and also only contributes to the drift term in the FFPE. The boson-boson interaction term is where for two position coordinates z and z # with the φ k (z) being any set of suitable orthonormal mode functions satisfying the condition dz φ k (z) * φ l (z) = δ k,l . A unitary change in these modes does not change δ C (z, z # ). The exact FFPE involves only first and third-order functional derivatives -there are no secondorder 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 ψ(z), ψ + (z) ∼ √ N , where N is the number of bosons present. As applying a further derivative gives a contribution of order 1/ √ N smaller than the previous term, discarding the third-order derivative terms is justified when N 1. The validity of the TWA is discussed in Refs. [21], [23], as well as in Section 2.2.4 and Appendix 7.

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]) and where the ψ(z, t), ψ + (z, t) are time-dependent stochastic field functions that now replace the original phase-space field functions ψ(z), ψ + (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 Γ a (t) = 0, Γ a (t 1 )Γ b (t 2 ) = δ a,b δ(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 ψ(z, 0) and ψ + (z, 0), as if the fields were like classical fields. However, ψ(z, t), ψ + (z, t) are non-classical and stochastic because the initial functions ψ(z, 0) and ψ + (z, 0) are members of an ensemble of stochastic fields, with a distribution chosen to represent the features of the initial quantum state ρ(0). This in turn leads to an ensemble of time-dependent stochastic fields ψ(z, t), ψ + (z, t) representing the time evolution of the quantum state, this being equivalent to the time evolution in the Wigner distribution functional W [ψ(z)].

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 If we choose ψ + (z, t) = ψ(z, t) * and ψ(z, t) = Φ c (z, t) the Ito SFE has almost the same form as the GPE. However, as well as Φ c (z, t) being non-stochastic, the quantity δ C (z, z) is absent from the nonlinear term in the GPE that allows for boson-boson interactions in the mean-field approximation. As dz Φ c (z, t) * Φ c (z, t) = N and dz δ 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 δ C (z, z) term in the Ito SFE will be relatively unimportant compared to the ψ + (z, t) ψ(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 ψ(z, 0) and ψ + (z, 0) should be more important than the δ C (z, z) term in treating quantum effects for interacting bosons. We define throughout a unity normalised condensate mode function ψ c (z, t) as

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.

Mean Energy
The expression for the mean energy H = T r( H ρ) (12) in terms of the stochastic field functions ψ(z, t), ψ + (z, t) is where and φ 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 nonstochastic. An expression for the mean energy is given below in Eq. (51) in terms of gravitational modes.

Position Probability Density
The number operator is given by which is the integral over an operator N (z) =Ψ † (z)Ψ(z) associated with the position probability density for the bosons in the system. We consider states which are eigenstates of N with eigenvalue N . The position probability density (PPD) is defined as the mean value of N (z) 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Ψ † (z)Ψ(z) in the , so that the position number density is then given either in terms of a functional integral or equivalently via a stochastic average = ψ(z, t) ψ + (z, t) − The stochastic average is the expression used in numerical calculations. Note the presence of the 1 2 δ C (z, z), but as the spatial integral of this term is of order 1 2 n M , whilst that of ψ(z, t) ψ + (z, t) is of order N c 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 ψ(z, 0) and ψ + (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.

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 where z, z # are two spatial positions.
, the QCF can also be expressed as a stochastic average via This QCF can be expressed in terms of natural orbitals χ 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) * . Then if we introduce an orthonormal set of mode functions φ k (z) we can define a matrix P with elements P k,l = . It is easy to show that P is Hermitian, so it can be diagonalised via a unitary matrix U such that P = U ∆U † , where ∆ 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 χ i (z, t) and their occupation numbers p i (see Eq.(21)) as where the natural orbitals are 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.

One-Body Projector
We can define a one-body projection operator (OBP) onto any single particle state |φ for an n body system in first quantisation as 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 1 bosons in the condensate mode the single-particle state would be taken as the initial condensate mode function with z|ψ c = ψ 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 M can be written in second quantisation form using an orthonormal set of mode functions φ k (z) and their annihilation, creation operators a k , a † k using standard procedures in which the field operators are expanded in terms of a set of orthonormal mode functions φ k (z) as Ψ(z) = k φ k (z) a k and Ψ(z) † = k φ k (z) * a † k . We have for the case of the condensate mode function where the completeness results 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 Ψ(z # ) † Ψ(z) times ψ c (z # , 0)ψ * c (z, 0). Hence, we see that the mean value for the one-body projector operator is which is the double space integral of ψ c (z # , 0)ψ * c (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 ψ(z, t) and ψ + (z, t) we get Note that the result involves the condensate mode function ψ 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 depletiondefined 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 ψ(z, t) and ψ + (z # , t) and the condensate mode function ψ c (z, 0) as a function of time. The term 1 2 δ C (z, z # ) would be small compared to the ψ(z, t) ψ + (z # , t) term for N c 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 assuming the initial condition of all bosons in mode ψ 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 since the initial condensate mode function ψ 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 where Φ 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) where 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 the mean-field position probability density.

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 ψ c (z, 0). One case of particular interest is where the mode function is one of the two Wannier-like modes Φ 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 M W i (t) and its mean value

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 ψ c (z, t) = Φ c (z, t)/ √ N obtained from the solution of the time-dependent GPE (10) and any set of orthogonal non-condensate modes ψ k (z, t) as where c 0 (t), c † 0 (t) are condensate mode annihilation, creation operators and γ 0 (t), γ + 0 (t) are the related stochastic amplitudes. Similar operators c k (t), c † k (t) and stochastic amplitudes γ k (t), γ + k (t) apply for the non-condensate modes.
The number of bosons left in the original condensate mode N C (t) is given by We can now derive an expression for the quantum depletion (QD) 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 By using Eq. (83) for the stochastic fields Ψ(z, t), Ψ + (z, t) the stochastic phase amplitudes for the condensate mode are given by Hence, using Eqs. (19) and (20) and 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 ψ c (z # , t)ψ c (z, t) * 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 Subsituting for the QCF in terms of stochastic fields gives has again been used.

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.

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 and are periodic φ k (z, t) = φ k (z, t + T ) as well as being orthonormal dz φ k (z, t) * φ l (z, t) = δ k,l at all times t. The Floquet frequencies ν 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 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 [ a k , The phase-space variables α k , α + k that represent the mode annihilation, creation operators are also time dependent -as are their related stochastic phase-space variables α k , α + 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.

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 wavefunctions 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 Expansion of the field function and stochastic field function in terms of these modes gives The phase-space variables η k , η + k that represent the mode annihilation, creation operators are time independent -whereas their related stochastic phase-space variables η k , η + 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 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 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 The derivation is set out in Appendix 8.

Evolution of Stochastic Phase Variables for Gravitational Modes
Coupled equations for the stochastic phase-space variables η k , η + k can be obtained allowing for the effect of the oscillating potential mg E z λ sin ωt and the boson-boson interactions. They can be written in the form 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 η k and η + k via the stochastic field functions, and their time dependence will ultimately depend on T and the 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
This requires a consideration of how to treat the preparation process, as will be explained below.

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 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 Φ (n+m)T (z, t + mT ) = Φ nT (z, t), which shows that for a given position z the Wannier mode Φ 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 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 φ 1 (z, t) and φ 2 (z, t). For effects involving a period sT we choose 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.

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 H of the energy, subject to constraints such as the mean boson number N 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 Φ c (z) obtained from the variational approach as the solution of a timeindependent GPE, but with the potential energy term now given by the time-independent trapping potential V trap (z) and including the chemical potential µ to allow for the constraint that N = 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 secondorder 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 Φ 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.

Condensate Wave Function, Condensate and Non-Condensate Fields
The condensate wave-function Φ c (z) is determined from the time-independent GPE where µ is the chemical potential. The condensate wave-function is normalised in terms of the number In Bogoliubov theory the field operatorΨ(z) is first written as the sum of the condensate field Ψ c (z) and a fluctuation (or non-condensate) field δ Ψ(z) aŝ where the condensate field and the fluctuation field are given by The condensate mode annihilation operator is c 0 and the condensate mode function is which is normalised to unity [32]. The equation for the condensate mode function is 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 ψ i (z), where the corresponding non-condensate mode annihilation operator is c i . The noncondensate modes are required to be orthogonal to ψ c (z), so we have the constraint It follows that the non-zero commutation result for the fluctuation field is 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 − 2 2m ∂ 2 However, in the present paper we will use Bogoliubov modes u k (z), v k (z) rather than Gross-Pitaevski modes.

Bogoliubov Modes
The fluctuation field operator is expanded in terms of annihilation, creation operators b k , b † k for Bogoliubov modes and associated mode functions The summation over k = 0 is to indicate that the condensate mode 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.
The commutation rules (67) for δ Ψ(z), δ Ψ(z # ) † together with the requirement that b k , b † k satisfy standard Bose commutation rules lead to the following conditions for the The inverse relations giving the b k , b † k in terms of the fluctuation fields are which is shown using Eqs (68), (70) and (66). Further conditions on u k (z), v k (z) follow using the commutation rules for b k , b † k and δ Ψ(z), δ Ψ(z # ) † and the orthogonality conditions (69) in conjunction with Eq. (71). These are 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 ω k we note that there are two mode functions involved.

Generalised Bogoliubov-de Gennes Equations
The Bogoliubov mode functions u k (z), v k (z) and the mode frequencies ω k are chosen to satisfy a generalised form [25] of the Bogoliubov-de Gennes (BDG) equations where k = 0. The differential operator L is defined by and the quantity C k is given by Note that the GPE (60) can be written L Φ c (z) = L ψ 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 Φ c (z), such as n C (z), Φ c (z) 2 and ψ c (z). Fortunately, the u k (z), v k (z) and the Bogoliubov frequencies ω 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 ω k are real, and that the u k (z), v k (z) satisfy the required orthogonality conditions (69) with the condensate mode function ψ c (z). The inclusion of the term involving C k is necessary to ensure orthogonality. By convention ω k is taken to be positive. The proofs make use of the Hermitian properties of L, The generalised BDG equations seem to have a solution for ω 0 = 0 with u 0 (z) = ψ c (z) and v 0 (z) = ψ * c (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 b k , b † k as the sum of uncoupled harmonic oscillators, one for each Bogoliubov mode. The form is given in Appendix 10 as Eq.(136).

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 u k , v k via the expressions Substituting into the generalised BDG equations we then obtain the standard BDG equations [21] ( where the orthogonality conditions (69) have been used. The quantity C k is now also given by in terms of these modes. A generalised orthogonality condition dz (ψ * c (z) u k (z) − ψ c (z) v k (z)) = 0 follows from the standard BDG equations, but this does not mean that the new modes u k , 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 u k , 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 / ω k . We find that

Field Operators -Condensate and Bogoliubov Modes
This gives the following expressions for the field operators after introducing the u k , v k into the expressions (61) for the fluctuation field. We havê 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 u k , v k involved are those obtained from the standard BDG equations (77).

Field Functions -Condensate and Bogoliubov Modes
As an alternative to Floquet or gravitational modes, the field functions Ψ(z), Ψ + (z) and the stochastic field functions can also be expanded in terms of Bogoliubov modes. The mode annihilation, creation operators b k , b † k are represented by phase-space variables β k , β + k , which are replaced by stochastic variables β k , β + k . Similarly, the condensate mode annihilation, creation operators c 0 , c † 0 are represented by phase-space variables γ 0 , γ + 0 , which are replaced by stochastic variables γ 0 , γ + 0 . Based on Eqs. (80) and re-introducing the original Bogoliubov modes we have 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.

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 ψ c (z, t) and any set of orthogonal non-condensate mode functions ψ i (z, t) -as in Eq. (63) for the fluctuation field operator. Thus, for the stochastic field functions we have (see Eq. (34)) The stochastic mode amplitudes are γ 0 (t), γ + 0 (t) for the condensate mode and γ i (t), γ + 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 γ 0 (0) = γ 0 , γ + 0 (0) = γ + 0 in view of ψ c (z, 0) = ψ c (z) -hence the same notation.

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 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 γ 0 , γ + 0 , β q , β + q for the Bogoliubov mode treatment of the BEC preparation can be translated into the choice of stochastic phase-space variables η k (0), η + k (0) for the gravitational mode treatment of the periodically driven BEC. Essentially, the η k (0), η + k (0) for the gravitational mode treatment are linearly dependent on the γ 0 , γ + 0 , β q , β + q .
We then see that we can write the initial stochastic fields for gravitational mode evolution as 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 (84) apply for the initial stochastic amplitudes α k (0), α + k (0) of Floquet modes with the replacement of the gravitational modes ξ k (z) by the Floquet modes φ k (z, 0), see Appendix 9.

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 β q (0), β + q (0) and γ 0 (0), γ + 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 . 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.

Density Operator for Initial State
The density operator is given by where This state is the product of vacuum states for the Bogoliubov modes and a Fock state with N bosons for the condensate mode.

First and Second-Order QCF for Initial State
We then find that where the terms involving v k (z) can be neglected since N 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 ψ c (z # ) * ψ 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 .

Wigner Distribution Functions for Initial State
For any given initial state, the stochastic averages of products of stochastic field functions ψ(z), ψ + (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Ψ(z),Ψ(z) † 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 γ 0 (0), γ + 0 (0), β q (0), β + 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 ρ C (0) for the condensate mode and a density operator ρ N C (0) for the non-condensate Bogoliubov modes, which are all in their vacuum states For this situation we can show that the Wigner distribution function W (γ 0 , γ + 0 , β, β + ) factorises into separate Wigner distribution functions W C (γ 0 , γ + 0 ) and W N C (β, β + ) for the condensate mode and the non-condensate Bogoliubov modes, where β ≡ {β 1 , ...β k , ..) and β + ≡ {β + 1 , ...β + 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 where the Wigner functions are normalised to unity -d 2 γ 0 d 2 γ + 0 W C (γ 0 , γ + 0 ) = 1 and d 2 β d 2 β + W N C (β, β + ) = 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, γ 0 β + k = γ 0 × β k . The result uses Morgan's approach [25] and involves first establishing a relation between the standard non-condensate mode operators c i , c † i that were introduced in Eq. (63) and the Bogoliubov mode operators b k , b † k which are given in matrix form (the convention differs from that in Ref [25]) as where the matrices U, V have elements given by This is a linear canonical transformation in which commutation rules are preserved. For the matrices, T denotes transverse, * denotes complex conjugation and † denotes the Hermitian adjoint. The bosonic commutation rules for both sets of non-condensate mode operators leads to the following matrix equations.
Further features of these matrices are set out in Appendix 13.

Stochastic Averages -Condensate Mode
For the condensate mode the density operator is given by (87) and as c 0 |N = √ N |N − 1 it is easy to show that the first and second-order QCF for normally ordered condensate mode operators are given by Since { c 0 c 0 } = c 0 c 0 , c † 0 c † 0 = c † 0 c † 0 and c † 0 c 0 = c † 0 c 0 + 1/2 we see that the stochastic averages are The distribution function generating the stochastic γ 0 and γ + 0 thus must produce a mean value of zero for the first order, and N + 1 2 for the non-zero second-order case γ + 0 γ 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 We then see that the non-zero stochastic average is This result is given in Ref [34]. If we choose γ 0 = Γ 0 exp(i φ 0 ) and γ + 0 = Γ 0 exp(−i φ 0 ), where Γ 0 is a normally distributed Gaussian random variable with a mean equal to N + 1 2 and a variance of 1/4, and we choose φ 0 to be a uniformly distributed random phase over 0 to 2π, then all the results in Eqs.(97), (98) and (99) are obtained. The Wigner distribution function for a Fock state is given in [34].

Stochastic Averages -Non-Condensate Modes
For the non-condensate mode the density operator is given by (90) and as b k |0 k = 0 × |0 k it is easy to show that the first and second-order QCF for normally ordered non-condensate mode operators are given by we see that the stochastic averages are The distribution function generating the stochastic β k and β + k thus must produce a mean value of zero for the first order, and 1 2 for the non-zero second order case β + k β k . It is of some interest to calculate the stochastic averages for the stochastic phase-space variables associated with the standard modes c i , c † i 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.

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.
In both the TWA and GPE numerics, the behaviour of these functions is deterministic and depends on the initial conditions ψ(z, 0) and Φ c (z, 0). The difference between the TWA and GPE is that in the former a stochastic ensemble of ψ(z, 0) is considered whereas in the latter just a single Φ c (z, t) is involved (see Sect.2.2). In the TWA numerics we find it convenient to choose stochastic field functions such that ψ + (z, t) = ψ(z, t) * and stochastic phase amplitudes such that η + k (t) = η k (t) * , 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 noncondensate 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.(84)).
The TWA numerical calculations are based on expanding the stochastic field functions in terms of time-independent gravitational modes, and during the calculations ψ(z, t) is known at each time point (see Sect. 2.3.2 and Eqs. (48), (52), (53)). At each stage ψ(z, t) and Φ 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 ψ(z, 0) and Φ c (z, 0) and then changing in shape and position as time evolves. If ψ(z, t) or Φ 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 ψ(z, t) or Φ 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 λ (in the oscillating frame). In carrying out these calculations we use (dimensionless) gravitational units: length l 0 = ( 2 /m 2 g E ) 1/3 , time t 0 = ( 2 /mg 2 E ) 1/3 , energy E 0 = mg E l 0 , where 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 where V is the potential in units of E 0 , and z, t are position and time in units of l 0 , t 0 respectively, with ω being the frequency in units of 1/t 0 . Henceforth the over-bars will be left understood.

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 λ and ω to be the same as theirs, namely λ = 0.12 and ω = 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 λ and ω result in two Floquet modes φ 1 and φ 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, ω 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 Φ 2 that is obtained from the two Floquet modes φ 1 and φ 2 . The Floquet frequencies are chosen to create a localised Wannier-like wave-packet which leads to (ν 2 − ν 1 ) being approximately ω/2. The various physical quantities for the simple case of a hard-wall mirror potential (V (z) = z for z ≥ 0 and ∞ 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 µ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, λl 0 /ω 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 .

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 Φ 1 and Φ 2 (Eq.  [1]. For|gN | less than for the critical interaction strength −0.006, the linear combination gives the Floquet mode φ 2 (z, t), which would have period T . When |gN | > 0.006, a more general linear combination of Φ 1 (z, t) and Φ 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 Φ 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 Φ 2 wave-packet is at the classical turning point z = h 0 . At times t = (k + 0.5)T (k = 0, 1, 2 . . .), the two wave-packets cross paths and interfere, again creating fringes.
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   (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 ω/2 (Fig. 3 (f)). Thus, discrete time-translation symmetry is broken at a critical interaction strength near gN ≈ −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 subharmonic peak in the FT at ω/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 ≈ −0.006, there is a small difference in the OBP (up to ≈ 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 ≈ −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.

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 Φ 1 , Φ 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 × 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 Φ 1 , Φ 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 Φ 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 ω/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 Φ 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 ≈ 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 ≈ −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 11 shows a logarithmic plot of the maximum quantum depletion due to quantum manybody 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 ≤ 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).

Quantum Depletion
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 (λ = 0), the TWA calculations indicate a relatively large quantum depletion, i.e., 40 − 220 atoms out of N = 600 for |gN | ≤ 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 ≈ 500 √ N , or at t/T ≈ 12, 000 mirror oscillations for N = 600. Figure 13 presents TWA calculations of the mean energy H 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.

Mean Energy of Atoms at Long Times
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.

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 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 meanfield 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 ≤ |gN | ≤ 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 ψ(z, t) and ψ + (z, t). For large N the δ C (z, z) term can be ignored, and the solutions for −g can be obtained from those for +g via the substitution ψ(z, t) = i ψ(z, t) # , ψ + (z, t) = i ψ + (z, t) # -which is possible for the double phase-space Wigner distribution W + where the field functions and stochastic field functions ψ(z, t) and ψ + (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.

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 nontrivial 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.

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). 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 H G = H − µ N , where H is the Hamiltonian, N is the number operator and µ is the chemical potential. We only consider number conserving states. As usual we consider a regime where the boson number N is large, N 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.

FFPE and Ito SFE -Extra Term
The density operator ρ satisfies the Liouville-von Neumann equation but for number conserving states [ N , ρ] = 0, so that we also have involving the grand canonical Hamitonian. The functional Fokker-Planck equation (see Eqs. (4), (5) and (6)) for the Wigner distribution functional W then has an additional term and this results in an extra term in the Ito stochastic field equations on the right hand side Consequently, from Eqs. (8) and (9), the full Ito stochastic field equations are now and ∂ ∂t ψ + (z, t)

Condensate and Non-Condensate Components
We first make the large N approximation of neglecting the terms δ C (z, z) in comparison to ψ + (z, t) ψ(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 Φ 0 (z, t) (or its conjugate) corresponding to the condensate mode and a stochastic term δ ψ(z, t) (or δ ψ + (z, t)) corresponding to the non-condensate modes. The stochastic terms δ ψ(z, t) or δ ψ + (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 If we substitute into Eqs. (109) and (110) and retain only the lowest order terms we find that This is the same as Eq. (2) in Ref. [1]. The stochastic non-condensate terms satisfy the equations and correct to the lowest order in δ ψ(z, t) and δ ψ + (z, t). Here, we see that the stochastic fluctuation fields δ ψ(z, t) and δ ψ + (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 δ ψ(z, t) and δ ψ + (z, t) are treated as stochastic quantities, whereas in timedependent Bogoliubov theory, the equivalent quantities would be time-dependent Bogoliubov mode functions u k (z, t), 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.

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 dz Φ * 0 (z, t)Φ 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 so minimising this with respect to changes in Φ 0 subject to the contraint leads to Eq. (112) above for Φ 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 Φ 0 . The advantage of this linkage is that if the condensate wave-function is expanded in terms of Wannier functions Φ i (z, t) which have periodicity sT 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ωt) that occur in manybody formulations of Floquet theory based on Shirley's paper [41], and which can be eliminated via a fast-rotating approximation.

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 timecrystal 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.

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Ψ(z), Ψ(z) † in each of the kinetic energy K, potential energy V and interaction energy 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 φ k (z) are also used.
We find that in the kinetic energy term where In the potential energy term we havê The interaction energy term is more difficult to deal with, but in this term we find that Hence, overall we have for the Hamiltonian 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 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 ψ(z), ψ + (z) then as stochastic averages involving the stochastic field functions ψ(z, t), ψ + (z, t) (see Sects 12.3.3 and 15.1.9 in Ref [22]). The final stochastic average expression for H is given in Eq (13).

Gravitational Modes
The mean energy expression (13) can be simplified via expanding the stochastic field functions ψ(z, t), ψ + (z, t) in some terms via gravitational modes. The potential energy is V (z, t) = mg E z − mg E λz cos ωt so the mean value of T + V is In the second and third lines we expand ψ(z, t), ψ + (z, t) using Eq.(48) and ∆δ C (z, z), δ C (z, z) via (7) and (14) using gravitational modes. This gives after spatial integration by parts where the defining equation (46) for the gravitational modes has been used, along with their orthogonality. Hence The mean value for the interaction energy U is left unchanged.

Position Probability Density and QCF -Floquet Modes
The position probability density in Eq. (18) can be expressed in terms of Floquet mode functions as 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 The quantity α k (t) α + l (t) − 1 2 δ k,l is the k, l element of a Hermitian matrix H, since α k (t) α + l (t) = (< a † l a k > + < a k a † l >)/2.

Evolution of Stochastic Phase Space Variables for Floquet Modes
By substituting the expressions in Eq (45) for ψ(z, t), ψ + (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) where the matrix L involves integrals of products of the Floquet mode functions, defined as 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 α k (0), α + 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 α k and α + 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 φ k (z, t) and the Floquet frequencies ν 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 where This alternative form is more convenient for numerical calculations. The non-linearity is embodied in the matrices D, D + . 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)).

Initial Conditions -Floquet Modes
These will be specified via the initial stochastic field which require knowing the Floquet mode functions at t = 0 and choosing a stochastic distribution of the α k (0), α + k (0) to match the quantum state that has been prepared in the trapping potential.

Grand Canonical Hamiltonian and Bogoliubov Approximation
The Hamiltonian 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 ρ is invariant under the U (1) symmetry group of phase changing unitary operators U (θ) = exp(−i N θ), it follows that [ N , ρ] = 0. Hence the evolution for ρ can be described by the grand canonical Hamiltonian K = H − µ N , with µ chosen so that N = N c .
In Bogoliubov theory the grand canonical Hamiltonian 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 c 0 with √ N c giving This is equivalent to writing the field operator asΨ(z) = Φ c (z) + δ Ψ(z), so the condensate field term is replaced by a non-operator field √ N c ψ c (z) = Φ 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 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) = Φ c (z) * Φ c (z). Note that a further approximation has been made -there are terms involving (δ Ψ(z) † ) 2 δ Ψ(z) 2 , (δ Ψ(z) † ) 2 δ Ψ(z), δ Ψ(z) † δ Ψ(z) 2 that have been discarded.
Since Φ 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.

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 In deriving Eq (136) the Hermitiancy properties of L, the reality of the ω k , the commutation rules for the b k , b † k , the biorthogonality conditions (72) and the orthogonality conditions (69) are all used, in addition to the generalised BDG equations (73).

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 |γ0 c for the condensate mode with amplitude γ 0 , and with each Bogoliubov mode k in its vacuum state |0 k . Thus For such a state we would have, if γ0 = √ N ĉ using c 0 |γ0 c = γ0 |γ0 c for Glauber coherent states. This state does obviously resemble a BEC with all N c bosons having the same wave function Φ c (z). However, it is not invariant under phase change transformations, so this is inconsistent with the requirement of ρ being phase invariant in order to treat evolution via the grand canonical Hamiltonian K. Nevertheless, a simple quantum state that is phase invariant can easily be constructed as a mixed state based on the |Φ γ0 . Writing γ0 = √ N c exp(iφ c ) we now consider the mixed state given by We then find thatΨ (z) |Φ γ0 = N c exp(iφ c ) ψ c (z) |Φ γ0 = exp(iφ c ) Φ c (z) |Φ γ0 Ψ (z) = 0 Ψ (z) †Ψ (z) = n C (z) = Φ c (z) * Φ c (z) 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.
The second order QCF results (149) also could be used to find expressions for stochastic averages involving pairs of Bogoliubov modes such as β + k (0) β m (0). However, this is more easily accomplished by considering the modes separately.

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 c i , c † i 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. (167) so that as β k = β + k = 0 it follows that γ i = γ + 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. (167) that where the stochastic averages of both sides have been taken. But from (103) we have so that γ + γ × γ γ + = 1 2 Evaluating these sub-matrices using (159) and (70) gives so that γ + i γ i = 1 2 + (V V † ) i,i , which is always greater than 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 γ i γ i = −(V U 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.

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.