Phase diagram and optimal control for n-tupling discrete time crystal

A remarkable consequence of spontaneously breaking the time translational symmetry in a system, is the emergence of time crystals. In periodically driven systems, discrete time crystals (DTC) can be realized which have a periodicity that is n times the driving period. However, all of the experimental observations have been performed for period-doubling and period-tripling discrete time crystals. Novel physics can arise by simulating many-body physics in the time domain, which would require a genuine realisation of the n-tupling DTC. A system of ultra-cold bosonic atoms bouncing resonantly on an oscillating mirror is one of the models that can realise large period DTC. The preparation of DTC demands control in creating the initial distribution of the ultra-cold bosonic atoms along with the mirror frequency. In this work, we demonstrate that such DTC is robust against perturbations to the initial distribution of atoms. We show how Bayesian methods can be used to enhance control in the preparation of the initial state as well as to efficiently calculate the phase diagram for such a model. Moreover, we examine the stability of DTCs by analyzing quantum many-body fluctuations and show that they do not reveal signatures of heating.


Introduction
Ever since the original conception of a time crystal in quantum many-body systems [1], there has been a growing interest to understand these objects theoretically [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17] as well as to realise them experimentally [18][19][20][21][22][23]. It turns out that in the model proposed in [1] such symmetry breaking cannot be observed if a system is initially in the ground state [24], while excited-state realizations are allowed [25]. Furthermore, time crystals cannot be observed for systems in the presence of long-range power-law interaction in their equilibrium state [26,27]. However, recent study shows [28] that the time crystal behaviour can also be observed in the ground state of a system with long-range interactions in the form of spin strings which are hard to realize in real experiments [29]. Popular realisations of time crystals involving excited states are DTCs which arise in periodically driven many-body systems. In such systems, the discrete temporal symmetry can be broken as a result of the inter-particle interactions. The first proposal of DTC involved a bouncing gas of ultra-cold atoms [2] and latter proposals involved spin-1/2 systems [3,4]. So far, the experimental realization of DTC has been performed in trapped ions [18] and nitrogen-vacancy centres in a diamond [19], where perioddoubling DTC and period-tripling DTC was observed respectively (see also [20][21][22][23]). Although there are few theoretical models to realize period n-tupling DTC [11,[30][31][32][33], where n > 3, there has yet been no experimental observation of them.
In this work, we study a model of ultra-cold bosonic atoms bouncing resonantly on an oscillating mirror [2] which has the potential of realizing n-tupling DTC, where n can be arbitrarily large. The motivation to physically realise n-tupling DTC lies in the fact that they provide a suitable platform to exhibit topological time crystals [34,35], temporal quasi-crystals [33,[36][37][38][39] as well as to demonstrate various nontrivial condensed-matter phenomena in the time domain [11,[40][41][42][43][44][45]. However, the experimental realization of this model involves multiple challenges as detailed in [11,46].
The formation of DTC in the system occurs due to sufficiently strong attractive interaction between atoms and it can be observed if an initial state is properly located with respect to the mirror position [11]. This requires a precise system control which cannot be guaranteed in each experimental realization. By evaluating the phase diagrams of this model, we examine the robustness of the DTC against perturbations of the initial atomic distribution. Moreover, we theoretically investigate the possibility of applying optimal control based on Bayesian optimization which could be performed on the experimental realisations of DTCs. We also address the question if the DTCs are stable against quantum many-body fluctuations which can result in heating of the system. Although the model we study allows arbitrarily large n, as a proof of principle, thorough analysis for period-doubling and period-quadrupling DTCs are addressed in this work. This will be hugely beneficial for real experiments which will most likely be carried out with a larger n to reduce possible atom losses.
The paper is organized as follows. In Section 2 we present the system of ultra-cold bosonic atoms bouncing on an oscillation atom mirror. The results concerning the phase diagrams and optimal control are shown in Section 3 and Section 4, respectively, while quantum many-body fluctuations of DTCs is presented in Section 5.

Theory: Model for period n-tupling discrete time crystal
We consider N ultra-cold bosonic atoms which are bouncing resonantly on an oscillating atom mirror. We assume that the system is strongly confined in the transverse directions and can be treated within the one-dimensional approximation. Atoms interact with each other via a contact potential whose strength g 0 is proportional to the s-wave scattering length and can be controlled by means of the Feshbach resonance mechanism [47]. We restrict the discussion to attractive interactions g 0 < 0. The system can be effectively described by a Bose-Hubbard model, Eq. (4), which is schematically presented in Fig. 1f. It turns out that for sufficiently strong attractive interactions, the system spontaneously breaks discrete time translation symmetry of the Hamiltonian and starts evolving with a period n-times longer than the period of the Hamiltonian (see Fig. 1a-e). Before describing the many-body phenomena in detail, it is useful to discuss the single-particle picture in order to better understand the concept of resonant driving in this model. (f) In the effective description, the system is described by Bose-Hubbard model, see (4).

Single particle problem
Consider a classical particle bouncing on an oscillating mirror with frequency ω and amplitude proportional to λ in the presence of a gravitational field. The Hamiltonian of the system in the frame moving with the mirror and in gravitational units has the form [2,11,[48][49][50] H 0 (z, p, t) = p 2 2 + z + λz cos ωt, z ≥ 0.
For λ = 0, the system is integrable and its motion is periodic with frequency Ω = ∂H 0 (I)/∂I, where H 0 (I) = (3πI) 2/3 /2 is the unperturbed part of the Hamiltonian (1) in the action-angle variables [11,[49][50][51]. In this picture, I is the classical analogue of the energy quantum number for an unperturbed particle. The distance of the mirror from the classical turning point is related to the frequency via the relation h 0 = π 2 /(2Ω 2 ).
In the presence of small mirror oscillations (λ 1), we are interested in the motion of a particle sufficiently close to a periodic orbit that is resonant with the driving, i.e., ω = nΩ, where n is an integer number. Because of the periodicity in time of the system, in the quantum description we can define the Floquet Hamiltonian H = H 0 −i∂ t which possesses T -periodic eigenstates (where T = 2π/ω) known as Floquet states [49,50]. To describe the motion of a particle close to the n : 1 resonant orbit one can apply the secular approximation [11,[49][50][51][52]. The resulting secular Hamiltonian indicates that in the frame moving along a resonant orbit, an atom behaves effectively like a particle in a periodic time-independent lattice potential with n sites and periodic boundary conditions and for n → ∞ a band structure in the quasi-energy spectrum emerges [11,40,[53][54][55][56]. We only consider the first quasi-energy band, therefore we construct n Wannier functions w i (z, t) localized in different sites of the periodic effective potential [11,40]. These Wannier functions, in the laboratory frame, are localized wavepackets w i (z, t) moving along a classical resonant orbit with period nT . Now we switch from single particle to ultra-cold atoms which fulfil an n : 1 resonance condition with the mirror motion.

Cloud of ultra-cold atoms
The many-body Floquet Hamiltonian of ultra-cold bosonic atoms which are bouncing resonantly on an oscillating atom mirror, in the Hilbert subspace of (nT )-periodic states, can be written in the form [2,11,40] where H 0 is the single-particle Hamiltonian (1) andΨ is the bosonic field operator. For a BoseEinstein condensate (BEC) all atoms occupy the same single-particle state and the many-body wave-function factorizes as φ 0 (z 1 , t)φ 0 (z 2 , t)...φ 0 (z N , t). In the mean-field approximation φ 0 (z, t) is a solution of the GrossPitaevskii (GP) equation [47] i∂ t φ 0 (z, t) = In order to get intuition about solutions of the GP equation that describe resonant motion of atoms let us restrict the analysis to the resonant single-particle Hilbert subspace spanned by the n localized Wannier wavepackets w i (z, t) of the first quasienergy band of the n : 1 resonantly driven system (1), see [2,11,40] for details. The mean-field energy functional corresponding to the many-body Floquet Hamiltonian (2), in the Wannier basis φ 0 (z, t) = n i=1 a i w i (z, t), has the form is the tunneling amplitude of atoms between neighboring Wannier wavepackets while U ij = 2g 0 nT nT 0 dt ∞ 0 dz |w i (z, t)| 2 |w j (z, t)| 2 for i = j and U ii = g 0 nT nT 0 dt ∞ 0 dz |w i (z, t)| 4 describe the strength of the effective interactions between atoms. This energy E is actually the mean-field quasi-energy per particle. Extrema of E are given by solutions of the GP equation (3) and can be found analytically for the case n = 2 [2], and numerically for n > 2. It turns out that if the strength of the attractive interactions is smaller than a certain critical value, i.e. |g 0 N | < |g cr N |, the mean-field solution corresponding to the minimal energy E is of the form φ 0 (z, [2]. However, when the interaction strength is larger than the critical value |g cr N |, φ 0 (z, t) is not a uniform superposition of w i (z, t), which means that the system chooses a periodic solution evolving with the period n times longer than the period expected from the symmetry of the Hamiltonian (see Fig. 1 for the case n = 4). The discrete time translation symmetry is broken and a period n-tupling time crystal phase forms. The solution for sufficiently large interaction is given by the single wavepacket φ 0 (z, t) ≈ w i (z, t). For this reason it is much better to realize the experiment in the regime |g 0 N | |g cr N | [11].

Challenges in the realization of a period n-tupling DTC
In this subsection we discuss the most important challenges in the realization of a period n-tupling DTC by means of ultra-cold atom system bouncing on the oscillating mirror.
In the laboratory it could be difficult to realize the hard-wall mirror that we have assumed in all theoretical analyses. However, if a realistic Gaussian shape mirror (that can be produced by a repulsive light-sheet) is used, the same time crystal phenomena can be realized as in the hard-wall case. [46]. In the experimental realization of a DTC, the initial distribution of BEC will be prepared as a harmonic oscillator ground state matching the Wannier state at a classical turning point [11]. However, the preparation of this initial atomic distribution is subject to experimental imperfections. For example, the atomic cloud released from the harmonic trap may have non-zero initial momentum, can be displaced from the classical turning point and may have non-Gaussian shape. The displacement of the cloud introduces detuning of bouncing gas of atoms from the resonant driving by the oscillating mirror. Moreover, there are typically shot-to-shot fluctuations in the position of the atoms in subsequent repetitions of an experiment. We expect that a sufficiently strong attractive interaction between atoms will compensate small displacements of the initial position of the wavepacket from the classical turning point. To show the robustness of DTC against perturbations of the initial atomic distribution, we determine the phase diagram as a function of displacement and interaction strength. However, if the displacement parameter is larger than the critical value and the initial position is not stable the mirror frequency needs to be corrected. For determining the optimal mirror frequency we have used the Bayesian optimization method.
It is also important to take into account the potential heating sources as well as possible atomic losses occurring in the system. Hence, we examine the stability of DTC on a long time scale against quantum many-body fluctuations which can result in heating of the system. In order to reduce atomic losses, it is better to choose a higher ratio of response period to driving period [11]. For a larger value of n the numerical simulations are time consuming. Therefore we also discuss the Bayesian optimization in the context of reducing numerical cost which would be essential for potential experiments, where a higher n is required [11,46].

Constructing the phase diagram for period-doubling and period-quadrupling discrete time crystals
We investigate the robustness of DTCs with regards to imperfect preparation of the initial state of the BEC. We start from an initial state which is the optimal Gaussian approximation of the mean-field solution of the DTC located exactly at the classical turning point h 0 above the mirror The parameterω 0 is the frequency of the harmonic trap where the BEC is initially prepared and it is chosen such that the width of the atomic cloud matches the width of the Wannier wavepacket w i (z, t) at the classical turning point [11]. We consider an n : 1 resonant driving of atoms by an oscillating atom mirror (see Sec. 2.1). Spontaneous breaking of the discrete time translation symmetry of the system, and consequent formation of the DTC, occurs for sufficiently strong attractive boson-boson interactions [2,11]. Then, the initially prepared wavepacket (5) is expected to be returning to the vicinity of the initial position after each period nT , as one can observe for the case n = 4 in Fig. 1.
It is convenient to introduce the quantum fidelity function where φ 0 (t) is the solution of the GP equation (3). Such a quantity can be recovered from an average particle-number distribution which can by measured as long as the overlap w 1 (t)|w 2 (t) ≈ 0 at the measurement time t. The fidelity (6) for n = 2 (thus for the period-doubling DTC) and for sufficiently strong attraction is presented in Fig. 2b. The Fourier transform of F (t) shows a single peak located at the half driving frequency ω/2 (Fig. 2d). This peak is related to the subharmonic response of the system which is the signature of the period-doubling time crystal [5]. However, for weak attractive interactions (which correspond to the symmetry preserving regime) after the tunneling time t = π/J T we observe transfer of atoms to the second wavepacket (Fig. 2a) which evolves along the same 2 : 1 resonant trajectory and is delayed (or advanced, depending on the point of view) by T with respect to the initial wavepacket. After another period π/J, atoms tunnel back to the initial wavepacket and this dynamics continues with  period 2π/J. Consequently, the Fourier transform of the fidelity function (6) reveals splitting of the Fourier peak around ω/2 (Fig. 2c). Such two separated peaks in the frequency domain are consistent with the beating in the plot of F (t). Note that the decay of the (2T )-periodic evolution due to tunneling of non-interacting atoms takes place even if there is no detuning from the resonant driving. It is in contrast to DTCs in spin systems where in the absence of any detuning from the perfect spin flip, time evolution of the spin systems is still periodic with period 2T [3][4][5][18][19][20][21].
As mentioned in Sec. 2.3 various experimental imperfections may cause fluctuations in the average momentum of atoms and of the frequency of the harmonic trap where the atomic cloud is initially prepared. Hence, we consider the initial states where the initial average momentum p 0 and the frequencyω are sampled from the uniform distributions in the intervals (−δp 0 , δp 0 ) and (ω 0 − δω 0 ,ω 0 + δω 0 ), respectively. Furthermore the initial location of the wavepacket (7) is taken to be h = h 0 + , where the displacement parameter denotes the detuning of the system from a perfect resonant driving. This displacement parameter and the interaction strength g 0 N constitute the space of parameters for which we determine the phase diagram of the DTC. We expect that even if = 0, sufficiently strong interactions are able to stabilize the evolution of the DTC.  To determine the critical value of for a given g 0 N , we analyze the fluctuations of the amplitude of the peak at ω/2 (for period-doubling DTC) and at ω/4 (for periodquadrupling DTC) in the Fourier transform of the fidelity (6) obtained in m random realizations of (7). We expect that the largest fluctuations of the peak amplitude can be observed near the critical point between the DTC regime and the symmetry unbroken regime [5]. Indeed, the variance of the peak amplitude shows a strong maximum at the critical point as clearly visible in Fig. 3a-b. Performing numerical simulations for different values of g 0 N we obtain the phase diagram depicted in Fig. 3c for the perioddoubling DTC. Similar phase diagram but for the period 4-tupling DTC is shown in Fig. 3d. In the latter case the maximum of the variance of the amplitude of the Fourier peak at ω/4 is used as the signature of the critical point. One can see that the phase diagrams are not symmetric with respect to = 0. It shows that the displacement with the case h > h 0 is favorable for the formation of a DTC, compared to the case h < h 0 . This is a consequence of the influence of the gravitational field. It should be noticed that the above approach is used in the regime where the symmetry breaking state is approximately given by a single wavepacket. Close to g cr N , the symmetry breaking states are superposition of two wavepackets with unequal weights. Therefore, in order to obtain g cr N , indicated by black dots in phase diagrams Fig. 3c-d, we have used the two-mode approximation (4).
In order to reduce the numerical burden, locating the critical , corresponding to a maximum of the variance, was turned into an optimization problem for which we have applied Bayesian optimization (see Appendix A). The positions of the variance peaks predicted at the end of the optimizations are marked by the dashed vertical lines in Fig.  3a-b. Using Bayesian optimization for 15 iterations to obtain the phase diagrams in Fig.  3c-d we have reduced by tenfold the computational time with efficiency around 70%.

Optimal control of the distance of the atomic cloud to the mirror
In an experimental realization of a DTC, the true distance h to the mirror may deviate more than the critical from h 0 (see Sec. 3) and will fluctuate around an average valueh in between each realization. Here we show that the mirror oscillations frequency could be adjusted, directly onto the experiment, to best account for this unknown average valueh.
For that purpose, we resort to Bayesian optimization in order to optimize the mirror frequency which allows one to control the true distance h. In the following we assume that the atomic cloud is prepared at a distance h randomly chosen, at each realization, from the uniform distribution in the range [h − δh,h + δh]. Because we do not know h we choose the frequencyω 0 of the harmonic trap optimal for h 0 . Furthermore,h is taken to be such thath − h 0 ≈ 0.3h 0 δh. In practice it is important to ensure thath is above h 0 because ifh − h 0 < −0.1h 0 , there is significant overlap of φ 0 (z, 0) with more than one Wannier state limiting the optimization procedure. To define the optimization problem one needs to specify a figure of merit to be maximized. In this case it is taken to be the fidelity function obtained after n periods of the mirror oscillations T = 2π/ω, where φ 0 (z, t = 0) is given in (7). At each step of the optimization this fidelity is averaged over 20 repetitions, and this average is fed into the optimizer. At the end of the optimization an optimal frequency ω opt is returned by the optimizer (see Appendix A). In Fig. 4, the average values of the fidelity are plotted as a function of discrete time t = k(nT opt ), where k is an integer number, T opt = 2π/ω opt and ω opt is the optimized frequency obtained after 15 iterations. It can be observed that the average fidelity is above 70% for a very long time for both the period-doubling (n = 2) and periodquadrupling (n = 4) cases. Without optimization the fidelity drops to almost zero after a few periods of the mirror oscillations. This shows the effectiveness of the Bayesian optimization technique based on a limited number of experimental repetitions to control the true distance h of the atomic cloud to the mirror.

Quantum many-body fluctuations of discrete time crystals
So far we have performed the analysis of the system within the mean-field approximation, i.e. according to the GP equation (3). These results allowed us to obtain the phase diagrams which determine how strongly one may perturb the system and still DTCs are stable. However, the mean-field approach assumes that time evolution of a Bose system is restricted to the many-body Hilbert subspace of product states, φ 0 (z 1 , t)φ 0 (z 2 , t) . . . φ 0 (z N , t). Interactions between bosons couple the product state subspace to the complementary space and can lead to decay of a DTC. We address this problem in the present section within the Bogoliubov approach [47] in the case of the period-doubling time crystal.
Before we switch to the Bogoliubov description, let us first comment on the behavior of ultra-cold atoms bouncing on an oscillating mirror, when we restrict to the manybody resonant Hilbert subspace [2,11]. In the case of the 2 : 1 resonant driving, the resonant Hilbert subspace is spanned by the two Wannier-like wavepackets w 1,2 (z, t), i.e. the bosonic field operator in (2) is restricted toΨ(z, t) ≈ w 1 (z, t)â 1 + w 2 (z, t)â 2 whereâ 1,2 are the standard bosonic annihilation operators. The many-body Floquet Hamiltonian (2) in the resonant subspace reads [2] with J and U ij similar to (4). The above HamiltonianĤ is identical to the Hamiltonian for bosons in a double well potential within the two-mode approximation [57] -in the present case the two modes are not time-independent functions but the (2T )-periodic Wannier wavepackets w 1,2 (z, t). It is known that the Hamiltonian (9) can be mapped to the Lipkin-Meshkov-Glick model [58] with a constant term omitted, where the spin operators read and γ = N (U 11 − 2U 12 )/J ∝ g 0 N . In the limit when N → ∞ but γ = const and γ < −1 (i.e. for sufficiently strong attractive interactions between bosons), the Lipkin-Meshkov-Glick model reveals a quantum phase transition where the Z 2 symmetry (the Z 2 symmetry means thatĤ commutes with e iπŜx ) of the Hamiltonian (10) is spontaneously broken. Then, all the eigenstates of the model for eigenenergies below the so-called symmetry broken edge −JN/2 reveal spontaneous breaking of the Z 2 symmetry. The spontaneous breaking of the Z 2 symmetry of the model corresponds to the spontaneous breaking of the time translation symmetry of ultra-cold atoms bouncing on an oscillating mirror which, in the basis of the (2T )-periodic Wannier functions w 1,2 (z, t), are described by the Floquet Hamiltonian (9) [2]. Thus, all Floquet many-body states of the Floquet Hamiltonian (9) with quasi-energies below the symmetry broken edge reveal perioddoubling time crystal behaviour. Even if in the resonant many-body Hilbert subspace the formation of the DTC is clear, there is still a question what happens in the full many-body Hilbert space of the system, i.e. when we take into account that interactions between bosons couple the resonant subspace with the complementary space? It can be addressed by applying the Bogoliubov approach. We use the particle-number-conserving version of the Bogoliubov theory [59] where the bosonic field operator is decomposed into the operator a 0 , corresponding to the condensate mode φ 0 (z, t), and the operator δΨ ⊥ (z, t) living in the orthogonal subspace, i.e.Ψ(z, t) = φ 0 (z, t)â 0 + δΨ ⊥ (z, t). The operator δΨ ⊥ (z, t) describes quantum many-body fluctuations around a many-body product state. Even within the Bogoliubov approach, it is not easy to calculate many-body Floquet states. However, it is much easier to investigate how a Bose system evolves in time if it is initially prepared as a perfect BEC. That is, when we start at t = 0 with all bosons in a many-body product state Φ(z, 0) = φ 0 (z 1 , 0)φ 0 (z 2 , 0) . . . φ 0 (z N , 0), we expect that interactions between bosons can lead to quantum depletion of a condensate. Initially we have a perfect condensate and consequently eigenvalues of the reduced single-particle density matrix, ρ(z, z ; t = 0) = Φ(0)|Ψ † (z, t = 0)Ψ(z , t = 0)|Φ(0) , are all zero except the one corresponding to a condensate mode φ 0 (z, t) which is equal to the total number of particles N . In the course of time evolution, bosons are being depleted from a condensate mode φ 0 (z, t) and start occupying other modes which is indicated by the fact that not only one eigenvalue of the reduced single-particle density matrix is nonzero. In the Bogoliubov approach the total number of bosons dN (t) depleted from the condensate is equal to the sum of the norms of v i (z, t) components of Bogoliubov modes In the particle-number-conserving version of the Bogoliubov theory [59], the modes evolve according to the following linear Bogoliubov-de Gennes equation where φ 0 (z, t) fulfills the GP equation (3),Q = 1 − |φ 0 (t) φ 0 (t)| andQ * = 1 − |φ * 0 (t) φ * 0 (t)|. Initial Bogoliubov modes allow us to define an initial many-body state of N bosons. If at t = 0 we choose [u i (z, 0), v i (z, 0)] = [χ i (z), 0] (where χ i |χ j = δ ij and χ i |φ 0 (0) = 0) and the Bogoliubov vacuum state as the initial state [60], we deal with the many-body state which is a perfect condensate [60], i.e. no bosons are initially depleted dN (0) = 0. The choice of χ i (z) is arbitrary but if we do our best and choose χ i (z) optimally adapted to a given problem, we will have to evolve a small number of the Bogliubov modes only in order to get the converged result for the total number of depleted bosons dN (t). In the case of the 2 : 1 resonant bouncing of ultra-cold atoms on an oscillating atom mirror, we start with φ 0 (z, 0) as the Gaussian state (5) and χ i (z) as eigenstates of the Hartree-Fock Hamiltonian H HF = − 1 2 ∂ 2 z + z + λz + 2g 0 N |φ 0 (z, 0)| 2 where contributions to the condensate mode are subtracted, i.e. χ i |φ 0 (0) = 0 and χ i (z) are corrected so that χ i |χ j = δ ij .
Evolving the condensate mode φ 0 (z, t) according to the GP equation (3) and the Bogoliubov modes [u i (z, t), v i (z, t)] according to the Bogolibov-de Gennes equations (13) we can obtain the reduced single-particle density matrix at any time t and diagonalize it, Figure 5a shows how the total depletion changes in time. It turns out that it is entirely determined by only one eigenmode of the single-particle density matrix, i.e. dN (t) = j dN j (t) ≈ dN 1 (t). After a long time of 1000 bounces of ultra-cold atoms on an oscillating mirror, i.e. at t f = 1999T , the total depletion is of the order of one atom only, dN (t f ) ≈ 1.6. It depends on the product g 0 N because the Bogoliubov-de Gennes equations (13) depend on g 0 N . In the experiment a typical total number of atoms is of the order N ≈ 10 4 which implies that for g 0 N = −0.02 chosen here, dN (t f )/N is negligible and the DTC is stable and resistant to quantum many-body fluctuations at least in a time scale that we have investigated here. The latter is much longer than expected duration of the experiment [11]. depleted from the condensate wavefunction φ 0 (z, t). Initially φ 0 (z, t) corresponds to the Gaussian wavepacket (5) and its time evolution φ 0 (z, t), according to the GP equation (3), describes the period-doubling time crystal. The total depletion of the condensate is dominated by a single eigenmode φ 1 (z, t) of the reduced single-particle density matrix (14), i.e. dN (t) ≈ dN 1 (t), because other eigenvalues dN j>1 (t) ≤ 6 × 10 −3 . Red dashed line shows the condensate depletion obtained with the help of the two-mode approach (9). Inset presents the total depletion dN (t) in the log-log-scale indicating the initial algebraic increase of the depletion. Right panel: probability densities of the condensate mode |φ 0 (z, t)| 2 (black solid line) and the dominant mode |φ 1 (z, t)| 2 (dotted-dashed line) at the final moment of the time evolution, i.e. at t = 1999T -both φ 0 (z, t) and φ 1 (z, t) are normalized to unity. The parameters of the mirror oscillations are: λ = 0.12 and ω = 1.4. The total evolution time 1999T corresponds to two tunneling periods of non-interacting atoms between the two Wannier wavepackets w 1,2 which is 2π/J where J = 7.26 × 10 −4 , cf. (9). The parameters of the two-mode Hamiltonian (9) are U/g 0 = 0.23 and U 12 /g 0 = 0.05, and the interaction strength is chosen so that g 0 N = −0.02. The two-mode results are obtained for N = 600 but they remain the same if N is greater. Figure 5b presents probability densities of the condensate mode |φ 0 (z, t f )| 2 and the dominant mode |φ 1 (z, t f )| 2 of the reduced single-particle density matrix (14). It turns out that | φ 0 (t f ± T )|φ 1 (t f ) | 2 ≈ 0.87 ± 0.02 what implies that atoms depleted from the condensate occupy the wavepacket that travels along the 2 : 1 resonant orbit but is delayed with respect to the condensate mode by the period T of the mirror oscillations. Thus, the many-body evolution of the system is restricted to two modes and these modes are similar to the modes w 1 (z, t) and w 2 (z, t) used to define the resonant many-body Hilbert subspace, cf. (9).
Let us compare the quantum many-body effects in the time crystal evolution obtained within the Bogoliubov approach and by means of the two-mode Hamiltonian (9). In the two-mode case, the perfect BEC with all atoms occupying the mode w 1 (z, t) corresponds to |N, 0 and it is the initial state we choose in the two-mode description. At t = 0 the quantum depletion of the condensate is zero but because the initial state is not an eigenstate of the Hamiltonian (9), the depletion increases in time which is shown in Fig. 5a. Despite the fact that in the Bogoliubov description, the initial condensate wavefunction φ 0 (z, 0) is not exactly the mode w 1 (z, 0) (i.e. φ 0 (z, 0) is the Gaussian approximation of the mode w 1 (z, t) only), the results for the quantum depletion obtained with the help of the both methods follow each other quite well. The results of the twomode approach correspond to N = 600 but they are the same for any N > 600 provided g 0 N = −0.02. The two-mode description allows us also to investigate what happens in an extremely long time scale. It turns out that at t ≈ 500T √ N , the depletion saturates at dN ≈ 0.02N and next, for much longer time evolution, the N -body system shows a revival and returns to the initial perfect BEC.
The initial BEC states used in the Bogoliubov description and in the two-mode approach are generic uncorrelated states of the resonantly driven many-body system. The presented results of the many-body time evolution of these states show that heating effects are negligible.

Conclusion
The intriguing possibility of studying novel quantum many-body phenomena in the time domain relies on realizing n-tupling DTC with large n, which so far has eluded any experimental observations. Adopting a model of ultra-cold atoms bouncing on an oscillating mirror that can be experimentally realized, we obtain the optimal conditions required for the initial state preparation of the system. Such optimization allows us to control and manipulate the system despite the experimental uncertainties and imperfections. The Bayesian method used here to obtain the phase diagram efficiently for small n can be naturally extended for large n-tupling DTC. This provides invaluable information for the experimentalists as it provides a well defined criterion for distinguishing the symmetry broken phase from the unbroken phase. Finally, our analysis of quantum many-body fluctuations that go beyond the meanfield approximation, clearly indicates that DTCs realized in our model do not show any signs of heating on long time scale -much longer than the duration of experiments. Thus, the integrability of the periodically driven single-particle system that reveals nonlinear resonances is inherited by the many-body counterpart if the interactions between particles are weak but still sufficiently strong to form DTCs. In summary, this work is definitely a step towards bridging the gap between theory and experiments on n-tupling DTC with the possibility of using optimal control in the experimental realisation of DTCs in ultra-cold systems.

Acknowledgements
We are grateful to Bryan Dalton, Peter Hannaford and Jia Wang for the fruitful discussion and valuable comments. In general, x can be a N -dimensional vector where N is the total number of parameters. In order to obtain the optimum solution x opt , one has to evaluate the figure of merit multiple times with F (x i ) representing the result of i th evaluation. The dependence of the figure of merit F (x) on x defines the optimization landscape which can often be non trivial. One way to perform this optimization task is by means of gradient methods. However, as it is the case here, analytical gradient of F (x) are not available, and approximations of these gradients by finite differences require extra numerical or experimental effort. Moreover this approach is limited by the presence of local extremum in the optimization landscape. Keeping this in mind, we resort to a non-gradient based optimization scheme, namely Bayesian optimization, which has been extremely efficient in terms of the number of evaluations needed to obtain close-to-optimum solutions [61][62][63][64] even in the presence of noise in the input parameters [65], or in the evaluation of the figure of merit [66]. Bayesian optimization is a technique that adopts a probabilistic approach towards optimization. It has essentially three steps: first, it approximates the unknown optimization landscape F (x) with well-behaved random functions f with prior distribution, that is before obtaining any evaluation, p(f ). Given The final step involves deciding, based on this model f , which parameters to evaluate next. A short-sighted strategy would consist on selecting the next set of parameters x i+1 where the model takes its maximal value. However, since the model f is only an approximation of the true optimization landscape F , there is also incentive to explore other regions of the parameter space where few evaluations have been recorded.
These two conflicting aspects, known as exploration-exploitation, are incorporated in Bayesian optimization by means of an acquisition function which values both the search for a maxima and also encourage exploration. The next set of parameters to evaluate is then chosen such that it maximizes this acquisition function. For example the upper confidence bound acquisition function, which was used in this work, is defined as µ f (x) and σ f (x) are respectively the mean and the standard deviation of the predictive distribution given in Eq. (A.2). The parameter k determines the exploration-exploitation balance with higher values of k corresponding to more exploration. More details about the specificity of building and updating the probabilistic models within Bayesian optimization can be found in Refs. [67][68][69][70].