Period-$n$ discrete time crystals and quasicrystals with ultracold bosons

We investigate the out-of-equilibrium properties of a system of interacting bosons in a ring lattice. We present a Floquet driving that induces clockwise (counterclockwise) circulation of the particles among the odd (even) sites of the ring which can be mapped to a fully connected model of clocks of two counter-rotating species. The clock-like motion of the particles is at the core of a period-$n$ discrete time crystal where $L=2n$ is the number of lattice sites. In presence of a"staircase-like"on-site potential, we report the emergence of a second characteristic timescale in addition to the period $n$-tupling. This new timescale depends on the microscopic parameters of the Hamiltonian and is incommensurate with the Floquet period, underpinning a dynamical phase we call 'time quasicrystal'. The rich dynamical phase diagram also features a thermal phase and an oscillatory phase, all of which we investigate and characterize. Our simple, yet rich model can be realized with state-of-the-art ultracold atoms experiments.

Introduction.-Symmetriespervade most fields of modern physics, ranging from nuclear physics to relativity and condensed matter physics.The spontaneous breaking of symmetries is the essential mechanism at the core of equilibrium phase transitions.Quite surprisingly, the possibility of systems breaking time-translational symmetry has been put forward only recently by Wilczek [1,2] who dubbed them 'time crystals' and triggered an extraordinary amount of excitement [3][4][5].Since time crystals have been shown to be impossible in quantum ground states [6], research has focussed on out-ofequilibrium conditions.Particularly successful has been the setting of Floquet systems, characterized by a time-periodic Hamiltonian H(t) = H(t + T ), for which the notion of a discrete time crystal (DTC) has been introduced in Refs.[7][8][9][10].A DTC is a system that, for some physical observable O with expectation value f (t) = ψ(t)|O|ψ(t) , in the thermodynamic limit and for a set of initial conditions |ψ(0) , features three properties [11]: (I) discrete time-translational symmetry breaking: f (t) = f (t + T ); (II) rigid subharmonic response: without fine-tuning f (t) shows oscillations with a period nT , i.e. period-n-tupling, that is the dynamics features a characteristic frequency 2π n with n integer ≥ 2; (III) persistence: the subharmonic oscillations extend up to infinite time.
In this letter we study the dynamical properties of bosons in a ring lattice.For a simple Floquet driving protocol of nearestneighbor hopping and local interaction, we find period-n DTCs (n ≥ 2), thermal, and oscillatory phases.Surprisingly, in the presence of a staircase-like on-site potential, a new dynamical phase emerges characterized by three incommensurate frequencies (Fig. 1).The dynamics is strictly aperiodic, but long-time ordered: a discrete time quasicrystal (DTQC) [27][28][29].
The idea of time quasi-crystallinity has recently appeared in the literature with various connotations and in various contexts such as dissipative classical systems [30], finite-size systems [31], and quasiperiodically driven systems [32][33][34].The term 'time quasicrystal' has also been adopted in an experiment with magnons in which one of the system's characteristic frequencies was incommensurate with the driving frequency [35] and in a model where a periodic repetition in time of a (possibly large, but finite) Fibonacci word was observed [25].
In striking contrast to the cases above, we consider a quantum, macroscopic, and periodically driven system and we find a DTQC phase that features a characteristic subharmonic frequency 2π n and two other, generally incommensurate, intrinsic frequencies ω 1,2  n .The first frequency is locked to a submultiple of the driving frequency, breaking discrete time-  3).In presence of interaction (U = 0.05) and a staircase-like potential (δ = 0.96) we observe a DTC (for = 0.01) and a DTQC (for = 0.05).The former is characterized by a single sharp peak locked at frequency 2π n , whereas the latter also features sharp peaks at two frequencies ω 1,2  n which are incommensurate with the driving frequency.
translational symmetry and robust to perturbations, similar to a standard DTC, whereas the latter two frequencies depend on the microscopic parameters.
The remainder of this letter is organized as follows.First, we introduce our model and study its solvable limits.We then propose a set of suitable diagnostic dynamical order parameters to characterize the system.We study the dynamics first in absence of an on-site potential and find a period-n DTC, a thermal, and an oscillatory phase.In the presence of a staircase-like on-site potential, we additionally find a DTQC phase.We demonstrate its rigidity, and find analytical expressions for the incommensurate frequencies 2π n and ω 1,2 n .Finally, we summarize our results and outline possible experimental implementations using state-of-the-art ultracold atom setups.
Model and integrable limits.-Weconsider a system of N bosons on a one-dimensional ring lattice with L = 2n sites governed by the Floquet Hamiltonian with period T = 1 [36] n jδ, that is a staircase-like on-site potential.The Floquet Hamiltonian (1) alternates hopping between odd nearest-neighbor links, hopping between even nearest-neighbor links, and onsite potential and two-body interaction.In the following t = 0, 1, 2, . . .denotes stroboscopic times.
We gain some intuition in the dynamical properties of the system by solving the Heisenberg equation of motion daj dt = i[H(t), a j (t)] ( = 1) in the two cases = 0 and U = 0.The case = 0 is at the core of a period-n DTC and is schematically illustrated in Fig. 1.Consider an initial product state in the site basis with n j bosons in an odd (even) site j.During the first part of the driving these particles move to site j + 1 (j − 1), during the second they move to site j + 2 (j − 2), and during the third the number of particles in each site is conserved (irrespective of U and {h j }).We thus find n j (t) = n j±2t (0) with + and − for even and odd sites j, respectively.The bosons in the ring lattice can be shown to be equivalent to a fully connected model of clocks with n-hands and of two species: clockwise and counterclockwise rotating (see SM).The clocks tick at every Floquet period, so that after a time t = n the system returns to its initial condition, a mechanism at the core of a period-n DTC.
In the non-interacting limit (U = 0), the Heisenberg equation of motion is linear and easily solved as a(t + 1) = F a(t), where a = (a 1 , a 2 . . ., a 2n ) T and where F is a 2n × 2n dimensional matrix.The dynamics is thus characterized by oscillations at the 2n frequencies corresponding to the phases of the 2n eigenvalues {λ j } of F .For instance, for = 0 we find doubly degenerate eigenvalues of the form λ 1,2 j = e iφ+i 2π n j (with φ some non-relevant phase), which correctly signals the period-n DTC with clock-like clockwise (counterclockwise) rotation of particles.
Dynamical order parameters.-Weare now interested in the dynamics away from the solvable limits = 0 and U = 0. To this end, we solve a semiclassical Gross-Pitaevskii equation of motion, which is obtained from the Heisenberg equation upon replacing the bosonic operators with c-numbers a → √ N ψ.In the thermodynamic limit of macroscopic occupation (N → ∞) for a fixed finite ring size L [3], considering a symmetry broken initial state with all the particles in site j = 1 (i.e.ψ j (0) = δ j,1 ), the semiclassical dynamics can either be chaotic, signaling quantum thermalization [37] or not, in which case we expect it to become exact [38].We indeed confirm this for our model explicitly using exact diagonalization and finite-size scaling (see SM).
It is now crucial to introduce suitable dynamical order parameters.To track the clock-like circulation of the particles we introduce a generalized imbalance on the odd sites and consider its Fourier transform A measure of the time crystallinity is given by whereas the chaoticness of the semiclassical dynamics is quantified by the distance d 2 between two slightly different initial states ψ(0) = (1, 0, 0, . . ., 0) T and ψ (0) = (1 − ∆, ∆, 0, . . ., 0) T with an arbitrary small ∆ = 10 −10 A growth of d 2 (t) to a finite value ∼ 1 corresponds to classical sensitivity to initial conditions and signals quantum thermalization [37].Finally, we consider the infinite time averages In the integrable limit = 0 we have a DTC, the phase of I grows linearly in time, Ĩ(ω) is peaked at ω = 2π n , Z(t) t = 1 and d 2 t ≈ 0. In the following, we will use these quantities to characterize the dynamical phases of the system.
Dynamical phases without on-site potential.-Westart by considering a vanishing on-site potential (δ = 0).In Fig. 2 we focus on the example L = 6, but similar results can be obtained for any L = 2n (see SM).We find that three dynamical phases are possible: (I) period-n DTC: for small ≈ 0 the clock-like circulation of the particles is rigidified by the interaction U = 0 and signaled by Z ≈ 1.The characteristic frequency of the imbalance I is locked to ω = 2π n , i.e. the frequency of the peak of Ĩ is robust against perturbations.This is the quintessential property of a DTCs: the interactions make the system subharmonic response rigid to mistakes in the driving, that is the range of corresponding to the DTC grows with U .(II) Oscillatory phase: for a large but small interaction U the system exhibits non-ergodic oscillations with a few characteristic frequencies signaled by sharp peaks in Ĩ(ω).In contrast to the case of a DTC, these frequencies are not locked, but rather depend on .(III) Thermal phase: for large and U classical chaos emerges, i.e. d 2 t ∼ 1.The three phases touch in the tricritical point = U = 0.
We emphasize that the observation of the DTC is not exclusive of the initial condition we considered here but is rather valid for a generic initial Fock state featuring a macroscopic imbalance, thanks to the underlying clock-like particle circulation.Moreover, in general our model does not require any disorder since in a fully connected system (to which our model can be mapped) disorder is not a necessary ingredient to prevent thermalization under Floquet driving [11,23].This is due to an extensive number of integrals of motion as explained in detail in the SM.
The realization of period-n DTC in a simple bosonic model in a ring lattice is the first major finding of this work.
Tilted lattice and discrete time quasicrystal.-Next,we show that tilting the lattice with a staircase-like potential (δ ≈ 1) favors a new DTQC phase.In presence of such a potential, we gain a clear intuition of the dynamics solving exactly the limit of U = 0 and δ = 1 (with generally = 0).In such a limit, the system's characteristic frequencies are linked to the eigenvalues of F and turn out to be (details in the SM) Only the first frequency is locked to a submultiple of the driving frequency, whereas the other two are generally incommensurate with it and are given by where p n ( ) is a trigonometric polynomial in with p n (0) = 1, e.g.p 2 ( ) = cos 2 and ω 1,2 2 = π ± 2 for L = 4.In the limit = 0 we find ω 1,2 n = 2π n , which consistently links back to the period-n DTC.In the thermodynamic limit -in our case corresponding to a macroscopic number of bosons N → ∞ but finite L -this regime represents a proper dynamical phase, which we will refer to as a DTQC [27,28].The discovery of a DTQC and of its characteristic frequencies (Eq.( 6)) represents the second major finding of this work.We emphasize that, differently from previous studies [32][33][34], the driving that we consider is perfectly periodic and characterized by a single frequency whereas the system response features a subharmonic frequency 2π n and two intrinsic frequencies ω 1,2  n that depend on the microscopic parameters of the system.
We show that these results are not an artifact of fine-tuning, but rather underline a proper dynamical phase.A non-zero interaction (U = 0) indeed rigidifies the DTQC.This means that, even if the frequencies ω 1,2  n depend on the model parameters (particularly on ), the characteristic subharmonic frequency 2π n does not shift.For δ ≈ 1 we observe the new DTQC phase, which is characterized by three main sharp  5), (c) the time crystal order parameter Z(t) in Eq. ( 4), and (d) the Fourier transform of the generalized imbalance Ĩ(ω) in Eq. (3).On-site interactions (U ) rigidify the DTC, with Z(t) ≈ 1 and Ĩ(ω) sharply peaked at ω = 2π/n for small (yellow lines).For larger and small interactions we find an oscillatory dynamical phase (OS), characterized by a few characteristic frequencies at which Ĩ(ω) is peaked (red lines).The thermal phase (TH) has Z oscillating chaotically, Ĩ(ω) with no dominant peak, and d 2 growing to a finite value ∼ 1 indicating sensitivity to the initial conditions (blue lines).The dynamical phases are identified as a function of and U via Z t (e) and d 2 t (f), which we compute as time averages over M = 10 4 Floquet periods.The dashed lines are the same in (e) and (f), the crosses correspond to the parameters considered in (b-d).We observe that interactions rigidify the DTC, which expands from = 0 for increasing U .peaks in Ĩ (Fig. 3a).One peak is locked to 2π n and does not change for and δ within a certain range, whereas the other two are well approximated by the non-interacting prediction of Eq. (6).We also still observe the DTC, oscillatory, and chaotic phases, which are reflected in Ĩ(ω) featuring a single peak locked to 2π n , few sharp peaks and no sharp peaks at all, respectively.In Fig. 3b we sketch the various phases as a function of δ and U for = 0.1 by looking at Z t = Ĩ( 2π n ) and d 2 t .The interaction U rigidifies the DTQC, which expands from δ = 1 for increasing U and is signaled by Z ≈ 0.5.For increasing interactions the system enters a thermal and/or a DTC phase.In Fig. 3c we characterize the system as a function of and U for δ = 0.9.We see that the interaction nec- essary to enter the DTC grows with .Eventually, if is too large, no DTC is possible.
Experimental implementations.-Anappealing feature of our proposal lies in its ultimate simplicity making it very natural for experimental implementation.Indeed, the Hamiltonian (1) only requires a number of basic ingredients such as nearest-neighbor hopping, local interaction and on-site potential, which are all readily available for ring-shaped optical lattices [39][40][41][42].The main difficulty is the discrete time dependent switch of tunnelling between alternating bonds.For small systems this is easily achievable if alternating bonds point in different real space directions.Alternatively, the even and odd site labels could be imprinted onto two different spin states.This would enable the discrete switch also on larger ring lattices.More concretely, the dynamics of the first and of the second third of the Floquet Hamiltonian would in this case correspond to a transverse field-induced spin-flip and to a intra-well barrier lowering, respectively.Furthermore, we have also checked that our results are robust against small undesired hoppings at odds with the alternating Floquet protocol.
Conclusions.-Inconclusion, we have studied the nonequilibrium properties of interacting bosons in a ring lattice.Considering a Floquet driving which alternates hopping on even and odd nearest-neighbor links, we induced a clock-like circulation of the particles whose direction depends on the parity of sites.Introducing a suitable set of dynamical order parameters and solving the exact dynamics in two inte-grable limits and a Gross-Pitaevskii equation across the whole parameter space, we identified several dynamical phases: a period-n DTC, an oscillatory phase, a thermal phase, and a new DTQC.The DTQC phase is characterized by a frequency which is locked to a submultiple of the Floquet frequency and by a second frequency which is incommensurate with it.Our work demonstrates a wide spectrum of non-equilibrium, nontrivial dynamical phases of bosons, which represent a unique opportunity for experimental investigation.
As natural for bosons in a finite-size lattice, the proposed system can effectively be mapped onto a fully-connected model.A natural question for future investigation regards then the fate of the various dynamical phases in systems breaking such a full-connectivity, e.g. in presence of long-range powerlaw interaction.Moreover, our work adds a totally new perspective to the multiple connotations of the still developing concept of time quasi-crystallinity, which certainly deserves further study.For instance, an intriguing question concerns the possibility for a system driven with two incommensurate frequencies ω 1 and ω 2 to respond at subharmonic frequencies ω 1 /n 1 and ω 2 /n 2 .
It is a pleasure to thank N. R. Cooper, J. P. Garrahan, M. Landini and A. Lazarides for useful discussions.A. P. acknowledges support from the Royal Society. A. N. holds a University Research Fellowship from the Royal Society and acknowledges additional support from the Winton Programme for the Physics of Sustainability.
Supplemental Material for "Period-n discrete time crystals and quasicrystals with ultracold bosons" Andrea Pizzi, Johannes Knolle and Andreas Nunnenkamp

I) DYNAMICAL EQUATIONS AND SOLVABLE LIMITS
We derive explicitly the quantum dynamical equation and solve it in the limit cases = 0 and U = 0.The Heisenberg dynamical equations for the annihilation operator a j at site j = 1, 2, . . ., L = 2n is obtained with a straightforward computation of the commutators of the Hamiltonian (1) with a j and reads ( = 1) where the upper and the lower signs are for odd and even j, respectively.In the following we implicitly assume a stroboscopic time t = 0, 1, 2, . . ., unless differently specified.Recalling J 1 = J 2 = π 2 + , we integrate the dynamics over the first two fractions of the Floquet driving to obtain whereas the third fraction of the Floquet driving can be integrated only in the two limit cases = 0 and U = 0.
For a fine-tuned hopping strength with = 0, we get and thus Moreover, since [H, n j ] = 0 for 2 < t < 3, we finally get where again the upper and lower sign refer to odd and even j, respectively.We thus obtain that n j (t = n) = n j (t = 0), that is the solvable limit = 0 is at the core of a period-n DTC.
In general, the system dynamics is characterized by 2n frequencies which are linked to the eigenvalues of F .Considering the staircase-like potential h 2j = h 2j−1 = 2π n jδ, these eigenvalues take a particular form in the limit δ = 1, which we are about to show and which is at the core of the DTQC phase.
Dividing both members of Eqs. ( 16) by ω j and introducing the Fourier transform that we rewrite in a compact form as where M k reads It is easy to show that det(M k ) = 1 for every k, so that M k admits an inverse M −1 k .Inverting Eq. ( 18), iterating it and exploiting periodicity in momentum space (ỹ which is itself a 2 × 2-dimensional eigenvalue problem for the eigenvalue λ n .The matrix M n M n−1 . . .M 1 has determinant 1, and its eigenvalues can therefore be written as −e ±i acos(− Tr 2 ) where Tr is its trace.For δ = 1, we thus finally get the eigenvalues of F to be where is in general a trigonometric polynomial in , of order 2n and such that p n (0) = 1.For instance, for L = 2n = 4 we obtain so that Tr = 2 sin( ) 2 − 2 cos( ) 2 , i.e. p 2 ( ) = cos(2 ).We plot p n ( ) for n = 2, 3, . . ., 10 in Fig. S1a, and report here its explicit expression for n = 2, 3, . . ., 6 As explained in the main text, the eigenvalues (21) are at the core of a DTQC dynamical phase.The characteristic frequencies ω 1,2 n are plotted in Fig. S1b.In the limit → 0 we get ω 1,2 n → 2π n , consistently recovering the DTC.

II) MAP TO A FULLY CONNECTED CLOCK MODEL
Our bosonic model is equivalent to a model of fully connected (i.e.infinite-range interacting) clock variables of two counterrotating species, as we show here in the limit = 0.For simplicity we adopt the inverse route: we present the clock model and map it back to the bosonic one, following in spirit the mapping that Ref. [23] carried out for clocks of a single species.refer to Appendix D of [23]).We say where P = i1,i2,...,i N Π i1,i2,...,i N is the symmetrization operator, with sum running over the possible permutations of 1, 2, . . ., N and with Π i1,i2,...,i N being the corresponding permutation operator.The bosonic operators are such that Since the Hamiltonian ( 33) is invariant under site permutations, it is easy to show that where n j = b † j b j is the number operator for the j-th bosonic mode.From Eq. ( 37), performing the sum over m, removing a constant term U and saying h 2j−1 = h o,j and h 2j = h e,j , we finally obtain As well, we find that the kickting operator acts as Summing up, we have therefore shown that the fully-connected clock model can be mapped into the model of bosons in a ring with 2n sites.The effect of the kicking operator is to rotate the bosons in the odd (even) sites of two steps in the clockwise (countercockwise) direction, whereas the clock Hamiltonian (33) maps into a local potential and two-body interaction for bosons.In the main, the kicking operator is then realized thanks to the sequential action of first and the second fractions of the ternary Floquet Hamiltonian.
The effective full-connectivity of the bosonic model has an important implication: the system is invariant under the permutation of clocks (1, 2), (2, 3), . . ., (N −1, N ), which corresponds to an extensive number N 2 of integrals of motion.Thanks to the existence of these integrals of motion, a priori the system does not necessitate any disorder to escape the fate of thermalization, making non-trivial dynamical phases such as DTCs possible.

III) ON THE VALIDITY OF MEAN FIELD
Throughout the main text, numerical results are obtained solving the Gross-Pitaevskii Equation (GPE), which is often referred to as the mean field (MF) limit of the Heisenberg equation.Here we discuss the regimes in which this is legitimate, and how to interpret the results when it is not.
In the considered thermodynamic limit (N → ∞ for a fixed L), the correct and well-established phase-space framework to deal with non-equilibrium dynamics is the Truncated Wigner Approximation (TWA) [38].This computes expectation values as averages over the classical GPE trajectories obtained for an ensemble of classical initial conditions, whose stochastic distribution is chosen consistently with the actual quantum initial condition.When the system is initialized in |ψ(0) = |N, 0, . . ., 0 , the initial condition for the GPE reads ψ j (0) = δ j,1 e iθj and, thanks to the gauge symmetry, can actually be considered without loss of generality to be ψ j (0) = δ j,1 .Since the initial condition no longer contains any degree of freedom, in this case the TWA is equivalent to a "single-shot" Gross-Pitaevskii Equation (SSGPE).This observation sounds very promising from a computational point of view, because it allows avoiding running the GPE multiple times, but it comes with a drawback: the SSGPE is not guaranteed to be accurate since the initial condition features non-macroscopically occupied sites.In particular, we expect the SSGPE to be inaccurate when chaotic.In the framework of TWA, chaotic semiclassical equations are in fact expected to underline thermalization [37], whereas the SSGPE displays in general persistent (in fact, chaotic) oscillations.We present results obtained within ED for L = 4 and in absence of on-site potential, showing evidence for the consistency of MF in the thermodynamic limit.We plot the time crystallinity order parameter Z (a1-a3), at stroboscopic times, for N = 5, 15, 25, comparing it with the MF result.We observe different trends corresponding to the different dynamical phases.In the DTC (a1) the agreement between ED and MF improves with for growing N and is rather good already for small N ; in the oscillatory phase (a2) the MF accuracy deteriorates in time, yet improves for increasing N : in the thermodynamic limit we expect the MF results to be exact at all times; in the thermal phase (a3) the MF is inaccurate irrespective of N , since thermalization quickly leads to ZED ≈ 0 whereas ZMF rather fluctuates chaotically.(b) We plot the measure of the error E of MF, see Eq. (40).For the DTC and the oscillatory phases, we find that E decays as a power of N (notice the logarithmic axes and the reference dashed lines ∝ 1/N, 1/N 1/2 ), whereas in the thermal phase E does not decay with N .

Exact diagonalization
To support the previous claims we compare the SSGPE with exact diagonalization (ED) computations for L = 4 sites, finite N ≤ 28 and vanishing on-site potential (δ = 0).For various values of and U corresponding to the DTC, oscillatory and thermal phases, in Fig. S2a1-a3 we plot Z(t).The DTC is characterized by Z ≈ 1 already for a relatively small N , whereas in the oscillatory phase the accuracy of the MF solution is worse and deteriorates in time.Nevertheless, in both cases, for increasing N the accuracy of MF improves, and in the thermodynamic limit N → ∞ (for a fixed, finite L) we expect MF to be exact up to arbitrarily large times.Conversely, in the thermal phase the MF solution is inaccurate irrespective of N .In this phase, Z M F fluctuates chaotically whereas Z ED saturates to a steady value, i.e. thermalizes, as expected [37].
To better interpret the results, we check the finite-size scaling considering the parameter which serves as a measure of the MF error with respect to ED for a given N .In Fig. S2b we show that E decays as a power of N in the DTC and oscillatory phases, whereas it does not decay in the thermal phase.This supports the idea for which, in the non-chaotic regime and in the thermodynamic limit, the SSGPE is exact.On the other hand, a failure of the SSGPE in the chaotic regime corresponds to the onset of quantum thermalization [37].

IV) LARGER RINGS
In the main paper we reported numerical results for L = 4, 6.However, our findings are valid for an arbitrary even number of sites L = 2n ≥ 4. To support this claim, here we corroborate our results studying a system with L = 8 sites, for which the time crystalline phases are characterized by period 4-tupling.To sketch the dynamical phases we look as usual at the height of the Fourier subharmonic peak Ĩ( 2π n ) = Z t and at the distance d 2 t between two initially very close copies of the system, which are shown in Fig. S3.The results that we find are completely analogue to the ones of the main text.

FIG. 1 .
FIG.1.Two entwined, counter-rotating clocks on a ring lattice.As a concrete example we show the n = 3 case.(a) Our Floquet driving scheme alternates (i) hopping between sites 2j and 2j − 1 (J1, single line), (ii) hopping between sites 2j and 2j + 1 (J2, double line), and (iii) on-site interaction and on-site potential, see Eq.(1).For J1 = J2 = π 2 (i.e.= 0) the particles in the odd (red) and even (blue) sites move clockwise and counterclockwise, respectively.(b) Fourier transform of the generalized imbalance Ĩ in Eq. (3).In presence of interaction (U = 0.05) and a staircase-like potential (δ = 0.96) we observe a DTC (for = 0.01) and a DTQC (for = 0.05).The former is characterized by a single sharp peak locked at frequency 2π n , whereas the latter also features sharp peaks at two frequencies ω1,2  n which are incommensurate with the driving frequency.

FIG. 2 .
FIG. 2.Period-n discrete time crystal, exemplified for n = L/2 = 3.For vanishing potential (hj = 0) we characterize the semiclassical dynamics via (b) the distance between two slightly different initial states d 2 (t) in Eq. (5), (c) the time crystal order parameter Z(t) in Eq. (4), and (d) the Fourier transform of the generalized imbalance Ĩ(ω) in Eq. (3).On-site interactions (U ) rigidify the DTC, with Z(t) ≈ 1 and Ĩ(ω) sharply peaked at ω = 2π/n for small (yellow lines).For larger and small interactions we find an oscillatory dynamical phase (OS), characterized by a few characteristic frequencies at which Ĩ(ω) is peaked (red lines).The thermal phase (TH) has Z oscillating chaotically, Ĩ(ω) with no dominant peak, and d 2 growing to a finite value ∼ 1 indicating sensitivity to the initial conditions (blue lines).The dynamical phases are identified as a function of and U via Z t (e) and d 2 t (f), which we compute as time averages over M = 10 4 Floquet periods.The dashed lines are the same in (e) and (f), the crosses correspond to the parameters considered in (b-d).We observe that interactions rigidify the DTC, which expands from = 0 for increasing U .

FIG. 3 .
FIG. 3. Discrete time quasicrystal.The on-site staircase-like potential leads to a DTQC, here exemplified for L = 4. (a) The dynamical phases are characterized by the Fourier transform of the generalized imbalance Ĩ(ω, ) for δ = 0.9, U = 2, and 500 Floquet periods.For a small | | we observe the DTC, corresponding to a single sharp peak locked at frequency 2π n .For larger we observe the DTQC, reflected by a sharp peak locked at frequency 2π n and two further peaks at frequencies ω 1,2 n ≈ π ± 2 (dashed blue lines).(b) The dynamical phases are characterized via the quantities Z t = Ĩ( 2π n ) and d 2 t (inset) computed over 10 4 Floquet periods and for = 0.1.The DTC, DTQC, and oscillatory phases correspond to Z ≈ 1, 0.5, 0, respectively, and thermalization is signaled by a finite d 2 t ∼ 1. Dashed lines serve as a reference and are the same in the main plot and in the inset.The DTQC is rigidified by interactions, expanding from δ = 1 for increasing U .For sufficiently large U , the system enters a DTC.(c) We characterize the system as a function of and U for δ = 0.9.The minimum interaction required to enter the DTC phase increases with | |, and no DTC is possible at all if | | is too large.
FIG. S2.Exact diagonalization results.We present results obtained within ED for L = 4 and in absence of on-site potential, showing evidence for the consistency of MF in the thermodynamic limit.We plot the time crystallinity order parameter Z (a1-a3), at stroboscopic times, for N = 5, 15, 25, comparing it with the MF result.We observe different trends corresponding to the different dynamical phases.In the DTC (a1) the agreement between ED and MF improves with for growing N and is rather good already for small N ; in the oscillatory phase (a2) the MF accuracy deteriorates in time, yet improves for increasing N : in the thermodynamic limit we expect the MF results to be exact at all times; in the thermal phase (a3) the MF is inaccurate irrespective of N , since thermalization quickly leads to ZED ≈ 0 whereas ZMF rather fluctuates chaotically.(b) We plot the measure of the error E of MF, see Eq.(40).For the DTC and the oscillatory phases, we find that E decays as a power of N (notice the logarithmic axes and the reference dashed lines ∝ 1/N, 1/N 1/2 ), whereas in the thermal phase E does not decay with N .

2 FIG
FIG. S3.Time crystalline phases in a L = 2n = 8-site ring.We characterize the DTC, DTQC, oscillatory and thermal phases with the parameter Z t (a,c) and distance d 2 t between two initialy very close copies of the system (b,d), considering up to 10 4 Floquet periods.