Discrete time crystals in Bose-Einstein Condensates and symmetry-breaking edge in a simple two-mode theory

Traditionally, spontaneous symmetry breaking (SSB) refers to a situation where the ground state or a thermal ensemble at a finite temperature of a many-body system is less symmetrical than its parent Hamiltonian. Breaking of different symmetries plays a profound role in many aspects of physics, including (spatial) crystal formation, magnetism, superconductivity and the origin of particle masses via the Higgs mechanism. However, spontaneous time-translational symmetry breaking had rarely been considered until Frank Wilczek proposed the controversial concept of a “time-crystal” (Wilczek, 2012). Wilczek’s original proposal, breaking of continuous time-translational symmetry in the ground state (or any thermal equilibrium state), was later rejected by the “no-go” theorem of Watanabe and Oshikawa (Watanabe and Oshikawa, 2015; Kozin and Kyriienko, 2019). Nevertheless, physicists have recently demonstrated that spontaneous discrete time-translational symmetry breaking (SDTTSB), i.e., discrete time crystals (DTCs), can exist in out-of-equilibrium systems, such as Floquet systems that are periodic in time with period $T$ (Sacha, 2015; Khemani et al., 2016; Else et al., 2016; Yao et al., 2017). Experimental evidence of time crystallinity has been reported recently in a variety of different platforms (Zhang et al., 2017; Choi et al., 2017; Rovny et al., 2018a; Pal et al., 2018; Rovny et al., 2018b; Smits et al., 2018; Liao et al., 2019; Smits et al., 2020). Discrete time crystals have developed rapidly recently and have attracted a lot of attention Else et al. (2017); Russomanno et al. (2017); Matus and Sacha (2019); Zhu et al. (2019); Sun et al. (2019); Machado et al. (2020); Natsheh et al. (2021); Pizzi et al. (2021). Reviews on the topic of time crystals can be found in Refs. (Sacha and Zakrzewski, 2018; Khemani et al., 2019; Else et al., 2020; Sacha, 2020).

However, a generic many-body system under periodic driving would normally keep absorbing energy and approach an infinite temperature state, a featureless state that cannot support SDTTSB (D’Alessio and Rigol, 2014; Lazarides et al., 2014; Ponte et al., 2015a). Therefore, the existence of a DTC relies on the prevention (or at least long-time suppression) of Floquet heating to stabilize the nonequilibrium quantum state. Indeed, several innovative mechanisms can help to avoid the heating problem in quantum many-body systems, including many-body localization (MBL) (Ponte et al., 2015b; Lazarides et al., 2015; Abanin et al., 2016), prethermalization (Kuwahara et al., 2016; Canovi et al., 2016; Machado et al., 2019) and, more recently, many-body quantum scars (Turner et al., 2018; Michailidis et al., 2020; Surace et al., 2021; Pizzi et al., 2020).

On the other hand, generalizing SSB to time-dependent Floquet systems where no well-defined ground state or any thermal equilibrium states exist requires a careful theoretical development (Watanabe and Oshikawa, 2015). It is now argued that a natural generalization of the equilibrium notion of SSB can be defined via the long-time steady state (Else et al., 2020). Following the nomenclature in Else et al. (2020), steady states in Floquet many-body systems are defined as those with expectation values of local observables that relax to constants at stroboscopic times $t=T,\ 2T,\ 3T,...$ In the case where time-translational symmetry is broken, we extend the definition to states with constant expectation values of local observables at $t=sT,\ 2sT,\ 3sT,...$ with $s$ being an integer. In a generic many-body system, any short-ranged correlated initial state would usually approach a steady state after some possible transient evolution. Therefore, in isolated quantum systems out of equilibrium, SSB can be defined as the situation where the steady state (in the thermalization limit) is less symmetrical than its parent Hamiltonian. In particular, a DTC refers to the case where the steady state has a period $sT$ that is a multiple $s$ of the drive period $T$. In the simplest case, for $s=2$, this amounts to a period-doubling or a subharmonic response at $0.5\omega\equiv\pi/T$.

A well-known property of Floquet systems is that a time-independent effective Hamiltonian, namely the Floquet Hamiltonian $H_{F}$, determines the stroboscopic ($t=T,\ 2T,\ 3T,...$) dynamics (Shirley, 1965; Sambe, 1973; Eckardt and Anisimovas, 2015). Therefore, the eigenstates and eigenvalues (namely Floquet states and Floquet energies, respectively) of $H_{F}$ or its equivalent Floquet evolution operator $U_{F}\equiv\exp(-iH_{F}T)$ also encode all the necessary information to define SDTTSB. Indeed, Refs. (von Keyserlingk et al., 2016; von Keyserlingk and Sondhi, 2016; Moessner and Sondhi, 2017) have shown that all eigenstates of $U_{F}$ for a MBL $s=2$ DTC come in pairs with eigenvalues that have a phase difference $\pi$, which is associated with the period-doubling of the steady state for systems starting with physically relevant initial states. In contrast, only one pair of eigenstates shows the $\pi$-pairing in a many-body quantum scar system studied in (Pizzi et al., 2020).

Here, we revisit the initial proposal Sacha (2015) that identified the possibility of realizing DTCs in a bouncing BEC under periodic driving, where not all but an extensive set of eigenstates with eigenvalues under a symmetry-breaking edge come in pairs. This initial work (Sacha, 2015) and some of the following studies applies a mean-field approximation (Giergiel et al., 2018, 2020) or a time-dependent Bogoliubov approximation (Kuroś et al., 2020), which might artificially preclude or underestimate Floquet heating in this system since only one or a few modes are included. Mean-field theories assume that there is no depletion from the condensate mode, and Bogoliubov theory assumes any depletion is small. As well as allowing for large depletion and allowing for quantum fluctuations, a multi-mode treatment is needed to examine the possibility of thermalization, as this could prevent DTC formation Khemani et al. (2019); Else et al. (2020). We recently investigated the case $s=2$ via a fully comprehensive multi-mode quantum treatment based on a phase-space many-body approach involving the truncated Wigner approximation (TWA) (Wang et al., 2021), which can include thermalization effects. However, thermalization is found to be absent in our system, and we find a robust sub-harmonic response for interactions stronger than some critical value $|g_{c}N|$ and lasting for a significant period of time, which is a practical criterion of DTC formation commonly adopted (Khemani et al., 2016; Else et al., 2016; Yao et al., 2017). Nevertheless, at present, the TWA is limited in the time regime that can be computed due to computation resource constraints.

Interestingly, we find for the parameter regime studied that only two (Wannier) modes were significantly occupied, suggesting that a many-body theory based on just two dominant modes should work well in this regime Wang et al. (2021). The two-mode model allows us to investigate the dynamics at much later evolution times. In this long-time regime, the system indeed relaxes to a steady state, where the dynamics are purified and the order parameter for time-symmetry breaking can be unambiguously determined. References (Sacha, 2015; Kuroś et al., 2020) also developed a beyond mean-field two-mode model. However, the approximation adopted in the derivation were not clearly justified. Here, we derive the two-mode model from a standard quantum field theory. We examine the two-mode approximation in comparison with the multi-mode TWA, and argue that the two-mode approximation for $s=2$ can be applied to a long evolution time. We also find that the spectrum of this model Hamiltonian determines the behavior and symmetry of the long-time steady state. We note that our two-mode model is not the same as two-mode treatments in mean-field theory (such as in Albiez et al. (2005)) where the condensate wave function is written as a linear combination of two mode functions. In contrast our two-mode model includes beyond mean-field effects and allows atoms to macroscopically occupy more than one mode. We found in our previous multi-mode TWA treatment that near the critical value of the interaction strength $|g_{c}N|$ (which corresponds to the onset of DTC formation), the mean-field theory does not provide a good description of the position probability density or one-body projector, nor can it account for depletion from the condensate mode (see Figs 6, 9(b), 9(f) and 11 in Ref. Wang et al. (2021)). The present two-mode model can account for all of these features.

The paper is organized as follows. In Sec. II we introduce the many-body model for a BEC of bosonic atoms bouncing resonantly on an oscillating mirror and investigate the validity of the two-mode model for the case of period-doubling. In Sec. III we present details of our two-mode model and calculations of the evolution of the system out to very long times. In Sec. IV we discuss the symmetry breaking in terms of a symmetry-breaking edge of the two-mode Hamiltonian. Our results are summarized in Sec. V. Details and some derivations of equations are given in the Appendices.

We consider $N$ bosons bouncing vertically on an oscillating mirror under strong confinement in the transverse directions, which can be regarded as a one-dimentional (1D) system. The quantum dynamics is determined by the many-body Schrödinger equation $i\hbar\partial t\left|\Theta(t)\right\rangle=\hat{H}\left|\Theta(t)\right\rangle$, where the Hamiltonian can be written

$$\hat{H}=\int dz\left[\hat{\Psi}(z)^{\dagger}H_{{\rm sp}}\hat{\Psi}(z)+\frac{g}{2}\hat{\Psi}(z)^{\dagger}\hat{\Psi}(z)^{\dagger}\hat{\Psi}(z)\hat{\Psi}(z)\right],$$

(1)

in terms of the field operators $\hat{\Psi}(z)$, $\hat{\Psi}(z)^{\dagger}$ for the annihilation, creation of a bosonic atom of mass $m$ at position $z$. Here, the single-particle Hamiltonian is given by $H_{{\rm sp}}=-\partial_{z}^{2}/2+V(z,t)$ using gravitational units for convenience, where the length, time and energy in gravitational units are $l_{G}=\left(\hbar^{2}/m^{2}g_{E}\right)^{1/3}$, $t_{G}=\left(\hbar/mg_{E}^{2}\right)^{1/3}$ and $E_{G}=mg_{E}l_{G}$ with $g_{E}$ being the gravitational acceleration. (Throughout this work, quantities shown in the figure are either dimensionless or are given in terms of gravitational units.) $g=2\omega_{\perp}a_{s}$ is the 1D coupling constant, where $a_{s}$ is the s-wave scattering length and $\omega_{\perp}$ is the oscillation frequency for the BEC atoms in a transverse trap. In the coordinate frame moving with the oscillating mirror (which can be transformed from the laboratory frame via a gauge transformation Buchleitner et al. (2002); Sacha (2015)), the temporally periodic potential is given by

$$V(z,t)=z(1-\lambda\cos\omega t),$$

(2)

with $z\geq 0$. Here, $\omega=2\pi/T$ is the driving frequency of the mirror, and $T$ is the corresponding period. The parameter $\lambda$ determines the driving amplitude. A formal method to solve the single-particle problem is the Floquet formalism, where the $T$-periodic Floquet eigenenergies $\epsilon_{\nu}$ and eigenstates $\phi_{\nu}(z,t)=\phi_{\nu}(z,t+T)$ can be defined as

$$\left[H_{{\rm sp}}-i\partial_{t}\right]\phi_{\nu}(z,t)=\epsilon_{\nu}\phi_{\nu}(z,t).$$

(3)

The single-particle classical motion under $H_{{\rm sp}}$ is chaotic for large $\lambda$, and only becomes regular with some suitable driving parameters and initial conditions. This regular motion can be recognized by the existence of regular resonance islands in the classical phase space that are located around periodic orbits with period $sT$, where $s$ is an integer Russomanno et al. (2017); Giergiel et al. (2018); Pizzi et al. (2021). Quantum mechanically, such regularity of the classical motion corresponds to the existence of $s$ special Floquet states $\phi_{\nu}(z,t)$, where $\nu=1,2,...,s$. Applying a unitary transformation to these Floquet states, one can construct $s$ Wannier-like states that are localized both in space and time and with temporal period $sT$ (Giergiel et al., 2018).

In this current work, we focus on the $s=2$ case with $\lambda=0.12$ and $\omega=1.4$ as an example (Kuroś et al., 2020). The Floquet states of interest $\phi_{1}(z,t)$ and $\phi_{2}(z,t)$ have quasi-energies $\epsilon_{1}\approx 0.410$ and $\epsilon_{2}\approx 1.109$. The two Wannier-like states are related to the two special Floquet states via

$$\begin{array}[]{l}\Phi_{1}(z,t)=\frac{1}{\sqrt{2}}\left[\phi_{1}(z,t)+e^{-i\pi t/T}\phi_{2}(z,t)\right]\\ \Phi_{2}(z,t)=\frac{1}{\sqrt{2}}\left[\phi_{1}(z,t)-e^{-i\pi t/T}\phi_{2}(z,t)\right]\end{array},$$

(4)

where one can verify that $\Phi_{\nu}(z,t)=\Phi_{\nu}(z,t+2T)$ and $\Phi_{1}(z,t+T)=\Phi_{2}(z,t)$. These Floquet states and Wannier-like states have been studied elsewhere (Kuroś et al., 2020; Wang et al., 2021) and are plotted in Fig. 1 of Ref. (Wang et al., 2021). The key property we note here is that at $t=0$, $\Phi_{2}(z,t=0)$ is a Gaussian-like wave-packet. Therefore, if we prepare a weakly interacting BEC confined in a harmonic trap with suitable initial position $\tilde{h}_{0}$ above the oscillating mirror and trap frequency $\tilde{\omega}_{0}$, almost all the atoms will initially occupy the mode $\Phi_{2}(z,t=0)$. In our previous TWA calculations in Ref. Wang et al. (2021), we chose $\tilde{h}_{0}=9.82$ and $\tilde{\omega}_{0}=0.68$, and found that at $t=0$, $N_{2}\approx 591$ atoms occupy mode $\Phi_{2}(z,t=0)$ out of a total atom number $N=600$. We used over 50 gravitational modes in our TWA treatment, and listed the other parameters in the Table 1 of Ref. Wang et al. (2021). This small difference between $N_{2}$ and $N$ is mainly due to the small mismatch between $\Phi_{2}(z,t=0)$ and a Gaussian wave-packet. We also find that the values of $\tilde{h}_{0}$ and $\tilde{\omega}_{0}$ do not need to be fine-tuned: a finite perturbation does not significantly reduce $N_{2}$ (Kuroś et al., 2020).

As pointed out in the Introduction, we have previously studied this system with finite interaction $g\neq 0$ using a multi-mode phase-space method, namely the truncated Wigner approximation (TWA) (Wang et al., 2021). Interestingly, we find that for $s=2$ during the whole evolution time $t\leq 2000T$ only these two Wannier-like modes are occupied. In Fig. 1 (a), we show the total occupation of these two Wannier-like modes $N_{s}(t)=N_{1}(t)+N_{2}(t)$, where $N_{i}\equiv N_{ii}$ is the occupation number of the mode $\Phi_{i}$. One can see that the fractional change of $N_{s}(t)$ as a function of time essentially has just small fluctuations around zero. There is also no obvious trend of it increasing during the time window investigated. (In stark contrast, many modes will be occupied if we switch off the driving for the same initial condition and interaction strength - see Ref. (Wang et al., 2021).) The fluctuations shown in this figure reflect the stochastic nature of the TWA calculations (and imply the actual value of $1-N_{s}(t)/N_{s}(0)$ might be smaller than the calculation uncertainty). In view of only two modes being important in the TWA calculations for $s=2$, in this work we apply a so-called two-mode approximation, where we project the Hamiltonian onto the Hilbert sub-space spanned by Fock states based on the two modes $\Phi_{1}$ and $\Phi_{2}$, which we name the stable-island Hilbert space. We emphasize here that, in contrast to fermionic systems that are limited by the Pauli exclusion principle, our bosonic modes can be occupied by an arbitrary number of atoms. Therefore, the two-mode model is a generic many-body problem for large particle number $N$, and the stable-island Hilbert space has a dimension of $N+1$.

While the two-mode model is described in detail in the next sections, here in Fig. 1 (b), we first show the excellent agreement of the TWA and the two-mode results for different interaction strengths $gN$. We emphasize that our TWA calculations Wang et al. (2021) show that a DTC forms for interaction strengths above $g_{c}N=-0.012$ for an initial state prepared in a harmonic trap with initial position $\tilde{h}_{0}$ and trap frequency $\tilde{\omega}_{0}$. Quantum fluctuations become significant near this critical interaction strength $g_{c}N$, and the TWA deviates strongly from the mean-field results. Nevertheless, the two-mode model is still in excellent agreement with the TWA in this case. The TWA in principle can include effects of many (much more than two) modes, and describe the many-body quantum dynamics exactly in the asymptotic $N\rightarrow\infty$ limit. On the other hand, although the two-mode model only includes two modes, it can describe the quantum evolution within the truncated stable-island Hilbert space exactly for any arbitary $N$. The excellent agreement between these two models gives us confidence that leakage of atoms to modes other than $\Phi_{1}$ and $\Phi_{2}$ is negligible during $t\leq 2000T$ and most of the atoms seem to be able to remain in the stable-island Hilbert space for a much longer time. The fact that most atoms are “trapped” in this stable-island Hilbert space for a long time breaks the ergodicity of the system and prevents Floquet heating. The underlying physics may be related to a quantum version of the Kolmogorov–Arnold–Moser (KAM) theory as pointed out in Ref. (Sacha, 2020), where the existence of a stable island (as tori in phase space) is robust against weak interactions (Lichtenberg and Lieberman, 1992). However, non-ergodicity in KAM theory can be destroyed by an Arnold diffusion after astronomically long time (Sacha, 2020), implying DTC in our system has a finite lifetime.

In the two-mode approximation, we assume the bosonic atom field operator can be expanded solely by the two Wannier-like modes

$$\hat{\Psi}(z)=\sum_{i=1,2}\hat{a}_{i}(t)\Phi_{i}(z,t),\ \hat{\Psi}^{\dagger}(z)=\sum_{i=1,2}\hat{a}_{i}^{\dagger}(t)\Phi_{i}^{*}(z,t),$$

(5)

where $\hat{a}_{i}(t)$ and $\hat{a}_{i}^{\dagger}(t)$ are the annihilation and creation operators of a boson in the time-dependent mode $\Phi_{i}(z,t)$. The time-dependence of the creation operators obeys $a_{i}^{\dagger}(t)=a_{i}^{\dagger}(t+2T)$ and $a_{1}^{\dagger}(t+T)=a_{2}^{\dagger}(t)$, which are determined by the properties of $\Phi_{i}(z,t)$, and similar rules apply to the annihilation operators. One can obtain a many-body basis set via the Fock state Dalton et al. (2015)

$$\left|n_{1},n_{2};t\right\rangle=\frac{\left[\hat{a}_{1}^{\dagger}(t)\right]^{n_{1}}}{\sqrt{n_{1}!}}\frac{\left[\hat{a}_{2}^{\dagger}(t)\right]^{n_{2}}}{\sqrt{n_{2}!}}\left|0,0\right\rangle,$$

(6)

where $n_{1}$ is the number of atoms in mode $\Phi_{1}$ and $n_{2}$ is the number of atoms in mode $\Phi_{2}$. Expanding the many-body state vector in this basis set gives

$$\left|\Theta(t)\right\rangle=\exp(-i\mu t)\sum_{n=0}^{N}b_{n}(t)\left|n,N-n;t\right\rangle,$$

(7)

where the total number of bosonic atoms $N$ is a good quantum number and $\mu$ is a suitable frequency (chosen below). These Fock states also satisfy $\left|n_{1},n_{2};t+2T\right\rangle=\left|n_{1},n_{2};t\right\rangle$ and $\left|n_{1},n_{2};t+T\right\rangle=\left|n_{2},n_{1};t\right\rangle$. We note that, while both the creation/annihilation operators and the Fock-state basis are time-dependent, the matrix elements of any direct product of creation and annihilation operators at the same time are time-independent, e.g., $\left\langle n_{1}+1,n_{2}-1;t\right|\hat{a}_{1}^{\dagger}(t)\hat{a}_{2}(t)\left|n_{1},n_{2};t\right\rangle=\sqrt{\left(n_{1}+1\right)n_{2}}$.

The time-evolution of the expansion coefficients $b_{n}(t)$ is determined by the many-body Schrödinger equation, in a vector form

$$\left[\underline{\widetilde{H}}(t)-i\hbar\partial_{t}\right]\vec{b}(t)=0,$$

(8)

where $\underline{\widetilde{H}}(t)$ is a time-dependent matrix with matrix elements $\widetilde{H}_{mn}(t)=\left\langle m,N-m;t\right|\hat{\mathcal{H}}(t)\left|n,N-n;t\right\rangle$. Here, the effective Hamiltonian operator is given by

$$\hat{\mathcal{H}}(t)=J\left(\hat{a}_{1}^{\dagger}\hat{a}_{2}+{\rm h.c.}\right)+\frac{1}{2}g\sum_{ijkl=1}^{2}U_{ijkl}(t)\hat{a}_{i}^{\dagger}\hat{a}_{j}^{\dagger}\hat{a}_{k}\hat{a}_{l},$$

(9)

where $J=(\mathcal{\epsilon}_{1}-\mathcal{\epsilon}_{2}+\hbar\omega/2)/2$ and $U_{ijkl}(t)=\int dz\Phi_{i}^{*}(z,t)\Phi_{j}^{*}(z,t)\Phi_{k}(z,t)\Phi_{l}(z,t)$. The phase factor $\mu$ occurring in the quantum state $\left|\Theta(t)\right\rangle$ is given by $\mu=N\epsilon/\hbar$, where $\epsilon=(\mathcal{\epsilon}_{1}+\mathcal{\epsilon}_{2}-\hbar\omega/2)/2.$ The derivation is given in Appendix A.1. With an understanding that the creation and annihilation operators will always act on the Fock basis and introducing time-independent matrix elements, we hereafter leave implicit the time-dependence of $\hat{a}_{i}$ and $\hat{a}_{i}^{\dagger}$. The only explicit time-dependence remains in the parameter $U_{ijkl}(t)$, which is periodic with period $\widetilde{T}=2T$ and $\widetilde{\omega}=\omega/2$. The effective Hamiltonian matrix $\underline{\widetilde{H}}(t)$ also shares the same periodicity, implying that Eq. (8) can be formally solved using a Floquet approach. A many-body Floquet state and Floquet energy can be defined as

$$\left[\underline{\widetilde{H}}(t)-i\hbar\partial_{t}\right]\vec{\mathcal{F}}_{\nu}(t)=\tilde{\mathcal{E}}_{\nu}\vec{\mathcal{F}}_{\nu}(t),$$

(10)

and the time evolution can be obtained via $\vec{b}(t)=\sum_{\nu}\mathcal{C}_{\nu}\vec{\mathcal{F}}_{\nu}(t)e^{-i\tilde{\mathcal{E}}_{\nu}t}$ and $\mathcal{C}_{\nu}=\vec{\mathcal{F}}_{\nu}^{\dagger}(t=0)\vec{b}(0)$ (see Appendix B). We emphasize here that this approach, namely the many-body Floquet (MBF) approach, is a full and exact many-body quantum calculation as long as only two modes are occupied and depletion to other modes is negligible.

Another time scale in the Hamiltonian is given by $t_{{\rm tunneling}}=1/J$ which describes the tunneling of a single particle between the two modes. $t_{{\rm tunneling}}$ typically is in the order of $500T$, much longer than $2T$, the period of $U_{ijkl}(t)$. An approximation to simplify and understand this problem is therefore to take a high-frequency expansion (HFE) of the Hamiltonian matrix $\underline{\widetilde{H}}(t)=\sum_{\delta}\underline{h}_{\delta}e^{-i\delta\tilde{\omega}t}$, where $\underline{h}_{\delta}=\int_{0}^{\tilde{T}}dt\underline{\widetilde{H}}(t)e^{i\delta\tilde{\omega}t}/\tilde{T}$, and keep only the lowest order $\underline{h}_{0}$ - which is a time-independent effective Hamiltonian matrix. The corresponding Hamiltonian operator $\hat{h}_{0}$ of $\underline{h}_{0}$ agrees with the one given in Ref. (Sacha, 2015), which can be expressed in the same form as Eq. (9) with $U_{ijkl}(t)$ replaced by the $2T$-average value $\bar{U}_{ijkl}=\int_{0}^{2T}U_{ijkl}(t)dt/2T$. However, we note that the time-independent Hamiltonian in Ref. (Sacha, 2015) is derived by expanding and truncating the Floquet Hamiltonian and field operators in an extended space-time Hilbert space. This expansion and truncation is thus an approximation that is effectively equivalent to our HFE approximation. Under the HFE approximation, the time evolution can be approximately given by $\vec{b}(t)=\sum_{\nu}c_{\nu}\vec{f}_{\nu}e^{-i\mathcal{E}_{\nu}t}$ and the initial condition $c_{\nu}=\vec{f}_{\nu}^{\dagger}\vec{b}(0)$. Here, $\mathcal{E}_{\nu}$ and $\vec{f}_{\nu}$ are eigenenergies and eigenstates of $\underline{h}_{0}$: $\underline{h}_{0}\vec{f}_{\nu}=\mathcal{E}_{\nu}\vec{f}_{\nu}$.

Discussions on SDTTSB usually refer to initial states without extensive many-body quantum entanglement, which are physically accessible by simple preparation schemes. As mentioned before, by preparing a BEC in a harmonic trap with suitable potential height and width, the initial condition can be well approximated by a product state $|0,N\rangle$. If we turn off the interaction in the preparation stage, the particles can tunnel from $\Phi_{2}$ to $\Phi_{1}$ with a tunneling rate $J$, and an additional relative phase can also be induced by some phase-printing scheme. Therefore, in principle, we can prepare a BEC initial product state, which in first quantization is

$$|\{\theta,\varphi\}\rangle=|\Lambda_{2}\rangle_{1}|\Lambda_{2}\rangle_{2}...|\Lambda_{2}\rangle_{N},$$

(11)

where all atoms occupy a single mode $\Lambda_{2}=\sin\theta e^{i\varphi}\Phi_{1}+\cos\theta\Phi_{2}$. Without loss of generality, we choose $\theta\in[0,\pi)$ and $\varphi\in[0,\pi)$.

To investigate the dynamical evolution, we study the observables $N_{ij}(t)=\langle\Theta(t)|\hat{a}_{i}^{\dagger}\hat{a}_{j}|\Theta(t)\rangle$, which play a similar role to the one-body reduced density matrix elements. For simplification of notation, we define $N_{i}(t)\equiv N_{ii}(t)$. Figure 2 shows a comparison between the MBF approach and HFE approximation for $gN=-0.006,\ -0.012$ and $-0.02$. The initial state is chosen as $|0,N\rangle$ wih $N=1000$. Figure 2 (a) (d) and (g) show perfect agreement between the HFE and MBF at relatively short time scales (0 to 2000$T$). At longer time scales $t>t_{{\rm relax}}$, Fig. 2 (b) (e) and (h) show that $N_{2}(t)$ relaxes to a constant $\langle N_{2}\rangle_{{\rm relax}}$ for the HFE. However, for an interaction strength near the critical value $gN=-0.012$, $N_{2}(t)$ has a relatively larger fluctuation. The MBF results here also show excellent agreement with the HFE, except for an almost negligible oscillation around the $\langle N_{2}\rangle_{{\rm relax}}$ that can only be seen in the insets for a large zoom-in scale. This small oscillation can be understood as micro-motion originating from the time-dependence of the many-body Floquet states in the MBF approach. After an even longer evolution time, Fig. 2 (c) and (i) show a strong quantum revival for both strong and weak interaction, but no revival can be seen for the critical case, Fig. 2 (f). We also note that for $gN=-0.02$, the quantum revival seems to be almost perfect, i.e. $N_{2}(t_{{\rm revival}})\approx N$, in contrast to the case for $gN=-0.006$ where $N_{2}(t_{{\rm revival}})<N$. Visible differences between the MBF and HFE show up in the quantum revival regime. Nevertheless, the HFE still accurately predicts the revival time $t_{{\rm revival}}$ and the overall revival profile. As one can see, the HFE is an excellent approximation. The HFE will also provide an insight into the symmetry breaking and give an analytical explanation of the time-evolution behavior, as shown in the next section.

While the quantum revival is interesting and indicates there is a relatively short transient period of deviation from an almost perfect temporal periodicity, the revival regime is much shorter than the steady-state regime. We discuss the quantum revivals with more detail in Appendix C and focus on the steady-state regime here. We emphasize that this relaxation value of $\langle N_{2}\rangle_{{\rm relax}}$ reflects the time-translational symmetry of the steady state. For example, the quantum correlation function (QCF) $P\left(z,z^{\#},t\right)={\rm Tr}\left[\hat{\Psi}^{\dagger}\left(z^{\#}\right)\hat{\Psi}(z)\left|\Theta(t)\right\rangle\left\langle\Theta(t)\right|\right]$ is given by

$$P(z,z^{\#},t)=\sum_{i,j=1}^{2}N_{ji}(t)\Phi_{i}(z,t)\Phi_{j}^{*}(z^{\#},t),$$

(12)

where $z$ and $z^{\#}$ refer to two different position coordinates. The position probability density at stroboscopic times $t=kT$ with $k=1,2,3...$ can be well approximated by $F(z,kT)=P(z,z,kT)\approx N_{1}(kT)\left|\Phi_{1}(z,kT)\right|^{2}+N_{2}(kT)\left|\Phi_{2}(z,kT)\right|^{2}$, where the cross-terms can be neglected since $\Phi_{1}$ and $\Phi_{2}$ are well separated localized wave-packets at $t=kT$ with $k$ being integer. We recall that $\Phi_{1}(z,kT)=\Phi_{2}(z,kT+T)$. Therefore, if $\langle N_{2}(t)\rangle\rightarrow\langle N_{2}\rangle_{{\rm relax}}=N/2$ as in Fig. 2 (b), $F(z,kT)=F(z,kT+T)$ has the same period as the Hamiltonian, which represents a symmetry-unbroken state. On the other hand, if $\langle N_{2}(t)\rangle\rightarrow\langle N_{2}\rangle_{{\rm relax}}\neq N/2$ as in Fig. 2 (h), $F(z,kT)\neq F(z,kT+T)$ and $F(z,kT)=F(z,kT+2T)$, indicating that the system’s state breaks the discrete time-translational symmetry, i.e., a time crystal is present. Figures 2 (h) and (i) predict that the time crystal survives for times out to at least 250,000 driving periods.

We can also define a one-body projection operator (OBP) $\hat{M}_{G}$ onto the Gaussian-like localized wave-packet $\Phi_{2}(z,0)$ as an observable (Wang et al., 2021). In the first quantization, the OBP is given by $\hat{M}_{G}=\sum_{i=1}^{N}(|\chi\rangle\langle\chi|)_{i}$, where $\langle z|\chi\rangle=\Phi_{2}(z,0)$ and $i$ lists individual particles as $i=1,2,...,N$. The expectation value of the OBP can be related to the QCF $\hat{\langle\Psi}\left(z^{\#}\right)^{\dagger}\hat{\Psi}(z)\rangle$ via $M_{G}(t)=\iint\mathrm{d}z\mathrm{~{}d}z^{\#}\Phi_{2}\left(z^{\#},0\right)\hat{\langle\Psi}\left(z^{\#}\right)^{\dagger}\hat{\Psi}(z)\rangle\Phi_{\mathrm{2}}^{*}(z,0)$. At stroboscopic times $t=kT$ with $k$ being integer, the OBP can be simplified as

$$M_{G}(kT)=\begin{cases}N_{1}(kT)=N-N_{2}(kT)&k=1,3,5,...\\ N_{2}(kT)&k=0,2,4,...\end{cases}$$

(13)

The physical picture is clear: $M_{G}(kT)$ measures how many atoms occupy the Gaussian-like wave-packet $\Phi_{2}(z,0)$ at stroboscopic times. When $k$ is even [odd], this Gaussian-like wave-packet coincides with $\Phi_{2}(z,kT)$ [$\Phi_{1}(z,kT)$]. Therefore, via the investigation of the behavior of $N_{2}(kT)$, we can obtain the evolution of the observable $M_{G}(kT)$, which reveals the temporal symmetry of the steady state. If $\langle N_{2}\rangle_{{\rm relax}}=N/2$, then $M_{G}(kT)\rightarrow N/2$ implies a $T$-periodicity. On the contrary, if $\langle N_{2}\rangle_{{\rm relax}}\neq N/2$, $M_{G}(kT)$ has $2T$-periodicity, indicating SDTTSB. One can define a frequency response observable via the Fourier transformation

$$m_{G}(f)=\frac{1}{K}\sum_{k=k_{0}}^{k_{0}+K-1}e^{-ifkT}\frac{M_{G}(kT)}{N},$$

(14)

where we usually choose $K=2048$ and $k_{0}=10^{4}$ so that $k_{0}T$ is in the steady-state regime. In particular, the sub-harmonic response is given by

$$|m_{G}(\omega/2)|\rightarrow\frac{|N-2\langle N_{2}\rangle_{{\rm relax}}|}{2N}=\frac{P_{G}}{2},$$

(15)

where $P_{G}=|\langle N_{2}\rangle_{{\rm relax}}-\langle N_{1}\rangle_{{\rm relax}}|/N$ is the population imbalance. $|m_{G}(\omega/2)|=0$ (or equivalently $P_{G}=0$) indicates the symmetry-unbroken phase, and hence $|m_{G}(\omega/2)|$ can be regarded as an order parameter.

Under the HFE approximation, the effective time-independent Hamiltonian satisfies the $\mathbb{Z}_{2}$ symmetry determined by the operator $\hat{P}_{12}=i^{N}\exp\left[i\pi(\hat{a}_{1}^{\dagger}\hat{a}_{2}+\hat{a}_{2}^{\dagger}\hat{a}_{1})/2\right]$ which interchanges the mode indices $1\leftrightarrow 2$. We find that $\hat{P}_{12}\hat{a}_{1}\hat{P}_{12}^{-1}=i\hat{a}_{2}$ and $\hat{P}_{12}\hat{a}_{2}\hat{P}_{12}^{-1}=i\hat{a}_{1}$. The corresponding effective Hamiltonian operator can be rewritten as

	$\displaystyle\hat{h}_{0}$	$\displaystyle=J\left(\hat{a}_{1}^{\dagger}\hat{a}_{2}+{\rm h.c.}\right)+gu_{T}\left(\hat{a}_{1}^{\dagger}\hat{N}\hat{a}_{2}+{\rm h.c.}\right)$		(16)
		$\displaystyle+\frac{1}{2}g\left[u_{I}\sum_{i=1}^{2}\hat{N}_{i}(\hat{N}_{i}-1)+4u_{N}\hat{N}_{1}\hat{N}_{2}\right]$
		$\displaystyle+\frac{1}{2}gu_{P}\left(\hat{a}_{1}^{\dagger}\hat{a}_{1}^{\dagger}\hat{a}_{2}\hat{a}_{2}+{\rm h.c.}\right),$

where $u_{T}=\bar{U}_{1112}$, $u_{I}=\bar{U}_{1111}$, $u_{N}=\bar{U}_{1212}$ and $u_{P}=\bar{U}_{1122}$ are the only four distinctive values of $\bar{U}_{ijkl}$ constrained by the $\mathbb{Z}_{2}$ symmetry (see Appendix A.2). The effective Hamiltonian $\hat{h}_{0}$ is invariant under the symmetry operator $\hat{P}_{12}.$ Here, the number operators are given by $\hat{N}_{1}=\hat{a}_{1}^{\dagger}\hat{a}_{1}$, $\hat{N}_{2}=\hat{a}_{2}^{\dagger}\hat{a}_{2}$ and $\hat{N}=\hat{N}_{1}+\hat{N}_{2}$. One might notice, under a mean-field approximation, that this effective Hamiltonian can be applied to investigate the bosonic self-trapping phenomenon, for example, in a BEC in a double-well potential Albiez et al. (2005). The connection of DTC formation and self-trapping was recognized in Ref. Sacha (2015). Here, we go beyond the mean-field approximation, and numerically investigate the whole eigen spectrum exactly. In our numerical investigation here, we have $J\approx 3.580\times 10^{-4}$, $u_{I}\approx 0.2237$, $u_{N}\approx 0.0519$, $u_{T}\approx-1.9\times 10^{-4}$ and $u_{P}\approx-4.3\times 10^{-6}$. Since the spectrum of $\hat{h}_{0}$ is bounded from both below and above, the maximum and minimum eigenenergy of the $\hat{h}_{0}$ can be obtained via the mean-field approach in the large $N$ limit. Neglecting the term associated with $u_{P}$ (which is much smaller than $u_{T}$, $u_{I}$ and $u_{N}$), this approach leads to

$$E_{{\rm max}}\approx E_{{\rm shift}}+|\tilde{J}|N/2$$

(17)

and

$$E_{{\rm min}}\approx\begin{cases}E_{{\rm shift}}-|\tilde{J}|N/2,&|gN|\leq|g_{b}N|\\ \frac{1}{2}gN^{2}u_{I}+\frac{\tilde{J}^{2}N}{gN(u_{I}-2u_{N})},&|gN|>|g_{b}N|\end{cases},$$

(18)

where $g_{b}N\equiv-2J/(u_{I}-2u_{N}+2u_{T})$. (We use the conditions $J>0$ and $gN<0$ here. We also focus on the case $J+gNu_{T}>0$ since $|gN|$ is small. See Appendix D for details.) Here, $E_{{\rm shift}}=gN\left(u_{I}N/2-u_{I}+u_{N}N\right)/2$ and its physical meaning will be clear below.

For a given total number $N$ (since $\hat{N}$ commutes with $\hat{h}_{0}$), the Hamiltonian can be mapped to the Lipkin–Meshkov–Glick (LMG) model (see Appendix A.3) as (Ribeiro et al., 2008; Mazza and Fabrizio, 2012; Kuroś et al., 2020)

$$\hat{h}_{0}\approx\tilde{J}\left(-\hat{S}_{x}+\frac{\gamma}{N}\hat{S}_{z}^{2}\right)+E_{{\rm shift}},$$

(19)

with $\hat{S}_{x}=(\hat{a}_{1}^{\dagger}\hat{a}_{2}+\hat{a}_{2}^{\dagger}\hat{a}_{1})/2$ and $\hat{S}_{z}=(\hat{a}_{2}^{\dagger}\hat{a}_{2}-\hat{a}_{1}^{\dagger}\hat{a}_{1})/2,$ which are bosonic spin operators. The LMG Hamiltonian can be applied to describe spin-systems with infinite-range interaction. However, we emphasize here that our system consists of bosonic atoms with zero-range contact interactions. When the two-mode and HFE approximations are applicable, the time-independent LMG Hamiltonian in Eq. (19) becomes an excellent approximation for the Floquet Hamiltonian of the underlying physical system, and determines the dynamical properties of the original system. The LMG Hamiltonian is invariant under the symmetry operator $\hat{P}_{12}.$ The parameters are given by $\tilde{J}=-2\left[J+gu_{T}(N-1)\right]$, $\gamma=gN\left(u_{I}-2u_{N}\right)/\tilde{J}$. In the case of large $N$ and finite $\gamma\propto gN$, a symmetry-broken edge exists for the LMG Hamiltonian Ribeiro et al. (2008); Mazza and Fabrizio (2012); Russomanno et al. (2017); Kuroś et al. (2020):

$$\mathcal{E}_{{\rm edge}}\approx-|\tilde{J}|N/2+E_{{\rm shift}}.$$

(20)

For $\mathcal{E}_{\nu}>\mathcal{E}_{{\rm edge}}$ the eigenenergies are non-degenerate whilst those for $\mathcal{E}_{\nu}<\mathcal{E}_{{\rm edge}}$ are essentially two-fold degenerate (see Fig. 3). From the expression for $E_{{\rm min}}$ and $\mathcal{E}_{{\rm edge}}$, one can observe that $E_{{\rm min}}=\mathcal{E}_{{\rm edge}}$ for weak interaction $|gN|<|g_{b}N|$, where $g_{b}N$ is defined above and Eq. (73). For strong interaction $|gN|>|g_{b}N|$, $E_{{\rm min}}<\mathcal{E}_{{\rm edge}}<E_{{\rm max}}$ and the eigenpairs below and above the edge have distinct features. We emphasise that, according to group representation theory, since $\mathbb{Z}_{2}$ is the symmetry group, the eigenstates for $\hat{h}_{0}$ would in general be non-degenerate and either even or odd under the symmetry operator, i.e., $\hat{P}_{12}|\nu,\pm\rangle=\pm|\nu,\pm\rangle$. However, for $|gN|>|g_{b}N|$, eigenstates that satisfy $\mathcal{E}_{\nu}<\mathcal{E}_{{\rm edge}}$ form near-degenerate pairs $|\nu,+\rangle$ and $|\nu,-\rangle$ of opposite symmetry, with an energy gap $|\mathcal{E}_{\nu}^{(+)}-\mathcal{E}_{\nu}^{(-)}|$ that is exponentially small in $N$. This energy gap is negligible for a large and finite $N$ in practice, and becomes exactly zero in the infinite $N$ (thermodynamic) limit. Therefore, we can choose a pair of symmetry-broken states that satisfy $\hat{P}_{12}\left|\nu\right\rangle_{U}=\left|\nu\right\rangle_{L}$ and $\hat{P}_{12}\left|\nu\right\rangle_{L}=\left|\nu\right\rangle_{U}$, which are given by $|\nu\rangle_{U}=(|\nu,+\rangle+|\nu,-\rangle)/\sqrt{2}$ and $|\nu\rangle_{L}=$$(|\nu,+\rangle-|\nu,-\rangle)/\sqrt{2}$. $|\nu\rangle_{U}$ and $\left|\nu\right\rangle_{L}$ have the same expectation energies $\mathcal{E}_{\nu}^{(U)}=\mathcal{E}_{\nu}^{(L)}$ and serve as a degenerate pair of symmetry-broken eigenstates. In contrast, any eigenstates with eigenenergies $\mathcal{E}_{\nu}>\mathcal{E}_{{\rm edge}}$ show no near-degeneracy and these states $\left|\nu\right\rangle_{S}$ are symmetry-unbroken. Figures 3 (a), (c) and (e) show the symmetry-breaking edge for $gN=-0.02$, $-0.012$ and $-0.006$, respectively. The vertical axis shows the minimum of adjacent gaps: $\min\left\{\Delta\mathcal{E}_{\nu}\right\}=\min(\mathcal{E_{\nu}}-\mathcal{\mathcal{E}}_{\nu-1},\mathcal{E}_{\nu+1}-\mathcal{\mathcal{E}}_{\nu})$, which is exactly zero if and only if there is a double-degeneracy. One can see that this quantity indeed becomes essentially zero at $\mathcal{E}_{{\rm edge}}$, which is denoted by the black vertical dash-dotted lines. The red dashed line on the right (left) indicates $\mathcal{E}_{{\rm min}}$ ($\mathcal{E}_{{\rm max}}$). For $|gN|\leq|g_{b}N|\approx 0.006$, only a symmetry-unbroken phase exists, since $\mathcal{E}_{{\rm edge}}\approx\mathcal{E}_{{\rm min}}$, and no near-degenerate pairs of energy levels occur.

A particularly useful observable to illustrate the importance of the symmetry-breaking edge is $\langle N_{2}\rangle_{\nu}$ and $\langle N_{1}\rangle_{\nu}=N-\langle N_{2}\rangle_{\nu}$, where $\langle O\rangle_{\nu}\equiv\left\langle\nu\right|\hat{O}\left|\nu\right\rangle$ and

$$\langle N_{2}\rangle_{\nu}=\begin{cases}\langle N_{2}\rangle_{\nu}^{(U)}\ {\rm or}\ \langle N_{2}\rangle_{\nu}^{(L)}&\mathcal{E}_{\nu}<\mathcal{E}_{{\rm edge}}\\ \langle N_{2}\rangle_{\nu}^{(S)}&\mathcal{E}_{\nu}\geq\mathcal{E}_{{\rm edge}}\end{cases}.$$

(21)

For symmetry-unbroken states, the permutation symmetry between modes ensures that $\langle N_{1}\rangle_{\nu}^{(S)}=\langle N_{2}\rangle_{\nu}^{(S)}=N/2$. On the other hand, for symmetry-broken states, we have $\langle N_{1}\rangle_{\nu}^{(U)}=\langle N_{2}\rangle_{\nu}^{(L)}$ and $\langle N_{1}\rangle_{\nu}^{(L)}=\langle N_{2}\rangle_{\nu}^{(U)}$, but in general $\langle N_{2}\rangle_{\nu}^{(U)}\neq\langle N_{2}\rangle_{\nu}^{(L)}\neq N/2$. Without loss of generality, we denote $\langle N_{2}\rangle_{\nu}^{(U)}>\langle N_{2}\rangle_{\nu}^{(L)}$ (hence the U- and L-branch). Figures 3 (b), (d) and (f) show $\langle N_{2}\rangle_{\nu}$ as a function of $\mathcal{E}_{\nu}$, with initial state $|0,N\rangle$ and $N=1000$. The black dash-dotted vertical line indicates the symmetry-broken edge, which illustrates the onset of bifurcation of $\langle N_{2}\rangle_{\nu}$ as a function of $\mathcal{E}_{\nu}$. We also show $\langle N_{2}\rangle_{{\rm relax}}$ by the purple pentagram symbol and the initial energy $E_{{\rm ini}}=\langle\Theta(0)|\hat{h}_{0}|\Theta(0)\rangle$ by the dashed vertical line. Numerically, $\langle N_{2}\rangle_{{\rm relax}}$ is calculated by averaging $\langle N_{2}\rangle_{{\rm relax}}=\sum_{k=k_{0}}^{k_{0}+K-1}\langle N_{2}(kT)\rangle/K$ over a time-window of $K=2048$ driving periods at $k_{0}T=10^{4}T$. We can see that the relaxation value lies on the $\langle N_{2}\rangle_{\nu}-\mathcal{E}_{\nu}$ curve of the corresponding branches, which can be understood via dephasing. In general, for an initial state $|\Theta(0)\rangle=\sum_{\nu}c_{\nu}|\nu\rangle$, the time-evolution of an observable is given by $\langle O\rangle(t)=\sum_{\mu\nu}c_{\mu}^{*}c_{\nu}\langle\nu|O|\mu\rangle e^{-i(\mathcal{E}_{\nu}-\mathcal{E}_{\mu})t}$. Typically, $c_{\nu}$ only has noticeable values around a small energy window close to the initial energy. If there is no degeneracy, after long enough time, dephasing leads to $\langle O\rangle_{{\rm relax}}=\sum_{\nu}p_{\nu}\langle O\rangle_{\nu}={\rm Tr}[\rho_{d}\hat{O}]$, where $\rho_{d}=\sum_{\nu}p_{\nu}|\nu\rangle\langle\nu|$ is sometimes called the diagonal ensemble, and $p_{\nu}=|c_{\nu}|^{2}$ is the projection probability of the initial state to the $\nu$’s eigenstate. Essentially, the contributions for $\mu\neq\nu$ cancel out for large $t$, as the phase factors in $\langle O\rangle(t)$ become more random. The initial energy can also be given by the diagonal ensemble $E_{{\rm ini}}=\sum_{\nu}p_{\nu}\mathcal{E_{\nu}}={\rm Tr}[\rho_{d}\hat{h}_{0}]$. When the initial energy is well above the broken edge, $\langle N_{2}\rangle_{\nu}=N/2$ gives $\langle N_{2}\rangle_{{\rm relax}}=N/2$ as shown in Fig. 4 (a). When the initial energy is below the broken edge, as shown in Fig. 4 (b) and (c), we find that for initial state $|0,N\rangle$, only the U-branch has noticeable projections, since the degeneracy is lifted in a single branch, the dephasing formula still works, and gives $\langle N_{2}\rangle_{{\rm relax}}=\sum_{\nu}p_{\nu}\langle N_{2}\rangle_{\nu}$.

To further explore the corresponding relation between the initial energy and the relaxation value, we investigate the case with an initial state $|\Theta(0)\rangle=|\{\theta,\varphi\}\rangle$ defined in Eq. (11) with all bosons in a mode given by $\Lambda_{2}=\sin\theta e^{i\varphi}\Phi_{1}+\cos\theta\Phi_{2},$ which has initial energy $E_{{\rm ini}}(\theta,\varphi)=\langle\{\theta,\varphi\}|\hat{h}_{0}|\{\theta,\varphi\}\rangle$. Since this initial state represents all atoms occupying the same single-particle state, the initial energy reduces to a mean-field expression

	$\displaystyle\frac{E_{{\rm ini}}(\theta,\varphi)}{N}\approx$	$\displaystyle J\sin 2\theta\cos(\varphi)+gN\left[\frac{u_{I}}{2}+\right.$		(22)
		$\displaystyle\left.u_{T}\sin 2\theta\cos(\varphi)+\frac{(2u_{N}-u_{I})}{4}\sin^{2}2\theta\right].$

Figure 5 shows $\langle N_{2}\rangle_{{\rm relax}}$ as a function of the initial energy $E_{{\rm ini}}(\theta,\varphi=0)$ for different initial states $|\{\theta,\varphi=0\}\rangle$, and $\langle N_{2}\rangle_{\nu}$ as a function of $\mathcal{E}_{\nu}$ for $gN=-0.1$. As one can see, in the deep regime in all branches (L, U and S), these two curves overlap, implying the dephasing mechanism leads to relaxation except very close to the critical point where finite-size effects become important. The inset shows the corresponding relationship between the initial state $\theta$ and the initial energy $E_{{\rm ini}}(\theta,\varphi=0)$. We also note that the initial states considered here are initial states with no multi-mode entanglement, as these are the initial states of interest for spontaneous symmetry-breaking physics.

In Fig. 6 (a), we investigate the subharmonic response $|m_{G}(\omega/2)|$ of the steady state as a function of $|gN|$ for different initial states $|\Theta(0)\rangle=|\{\theta,\varphi=0\}\rangle$, which changes from $0$ to a finite value abruptly around a critical $|g_{c}N|$ indicated by the thin vertical lines. These critical values for different $\theta$ can be obtained by equating $E_{{\rm ini}}(\theta,\varphi=0)$ and $\mathcal{E}_{{\rm edge}}$ as shown in Fig. 7 (a). The inset of Fig. 6 (a) shows the frequency response $|m_{G}(f)|$ for $f$ close to $\omega/2$ and the initial state $|\Theta(0)\rangle=|\{\theta=0,\varphi=0\}\rangle=|0,N\rangle$. We can see that a single sharp peak appears abruptly when $|gN|\geq|g_{c}N|$. However, in realistic experiments, one might not be able to access such very long evolution times. In Fig. 6 (b), we show that the transient behavior at relatively short times can be regarded as a precurser. One can see that, while not as smooth and clean as Fig. 6 (a), $|m_{G}(\omega/2)|$ still shows an abrupt change to a much larger value near $|g_{c}N|$. In the inset, one can see that the frequency response of the transient states for $\theta=0$ at short times shows a double-peak structure for $|gN|<|g_{c}N|$ and which changes to a single-peak structure for $|gN|\geq|g_{c}N|$, the same as the observation in (Wang et al., 2021). In the steady-state regime, we now have the intriguing situation where, as the magnitude of the interaction is increased from that just below the critical interaction $|g_{c}N|$ to just above $|g_{c}N|$, the period of the bouncing atom cloud changes from period $T$ to period $2T$, to form a DTC. Although the numerical results presented here are calculated for a finite $N=1000$, we believe the conclusions remain valid in the thermodynamic limit ($N\rightarrow\infty$ and finite $gN$). An analysis of the finite size effect is given in Appendix E.