Condensed Matter Physics in Time Crystals

Time crystals are physical systems whose time translation symmetry is spontaneously broken. Although the spontaneous breaking of continuous time-translation symmetry in static systems is proved impossible for the equilibrium state, the discrete time-translation symmetry in periodically driven (Floquet) systems is allowed to be spontaneously broken, resulting in the so-called Floquet or discrete time crystals. While most works so far searching for time crystals focus on the symmetry breaking process and the possible stabilising mechanisms, the many-body physics from the interplay of symmetry-broken states, which we call the condensed matter physics in time crystals, is not fully explored yet. This review aims to summarise the very preliminary results in this new research field with an analogous structure of condensed matter theory in solids. The whole theory is built on a hidden symmetry in time crystals, i.e., the phase space lattice symmetry, which allows us to develop the band theory, topology and strongly correlated models in phase space lattice. In the end, we outline the possible topics and directions for the future research.


Wilczek's time crystal
The idea of time crystals was proposed by Frank Wilczek in 2012 [1,2]. He raised the question that whether continuous time-translation symmetry (CTTS) might be spontaneously broken in a closed time-independent quantum mechanical system, namely, some physical observable of the quantum-mechanical ground state exhibits persistent periodic oscillation. Wilczek also noticed that the Heisenberg equation of motion for an operator, i.e., for any eigenstate Ψ E of H, seems to make the possibility of quantum time crystals implausible. Nevertheless, he proposed a model consisting of particles with contact attractive interaction on a ring threaded by an Aharonov-Bohm (AB) flux, and found that this model can have the "ground state" of a rotating soliton, which satisfies the definition of time crystals. However, the spontaneous breaking of CTTS is against the energy conservation law. The rotating soliton generates a time-dependent force field on the ring. Imagine that we put a particle at some point on the ring when the force exerted by the soliton is weak. As the soliton is approaching the point, the particle will be accelerated and absorb some energy from the rotating soliton. In this way, some positive net energy is produced without lowering the energy of the rotating soliton as the ground state has the lowest energy by definition. Wilczek's ring is actually a new version of perpetual motion machine, which implies there must be something wrong in his analysis. In the history of science, although all the proposals of perpetual motion machines failed, several famous ones, e.g., Maxwell's demon [3] and Brownian ratchet [4], played significant roles in advancing science and still inspire new work today. Disproving Wilczek's idea needs serious and rigorous mathematical proof, which may deepen our understanding of physics laws.

No-go theorem
Bruno first pointed out Wilczek's rotating solution is not the true ground state for the model he proposed [5]. The correct ground state is actually a stationary state without breaking CTTS [6]. However, it is too early to disprove Wilczek's general idea of time crystals by precluding his soliton model. In the original work [1], Wilczek started his analysis from considering a superconducting ring threaded by a magnetic flux that is a fraction of the flux quantum. If Cooper pairs are free of interaction, the ground state is indeed a constant uniform superconducting current along the ring. But if the Cooper pairs have some kind of interaction (not only restricted to the contact form), one can ask whether a charge density wave (or Wigner crystal) can appear and rotate as the ground state. Actually, this rotating change density wave was soon proposed with trapped ions after Wilczek's original work [7]. However, a detailed analysis showed that the rotating charge density wave cannot be the ground state [8]. A "no-go theorem" was proved later by Bruno [9], which rigorously rules out the possibility of spontaneous ground-state (or thermal-equilibrium) rotation for a broad class of systems.
Spontaneous symmetry breaking occurs in the thermodynamic limit, where the system size N (e.g., the number of particles in Wilczek's model or the number of spins in Ising or Heisenberg models) goes to infinity. A general method to treat the thermodynamic limit correctly was given by Bogoliubov [10]: first calculating physical quantities for the system with finite size N and perturbative symmetry breaking potential v, then taking the limit N → ∞ before taking the limit v → 0. Following Bogoliubov's prescription, Bruno found that, in the presence of an external symmetrybreaking potential rotating at angular velocity ω, the nth-level energy of the system in the static frame is where φ is the flux threaded through Wilczek's ring and I φ,n is the moment of inertia for the nth energy level of the system. Bruno's crucial finding is that I φ,0 > 0 for the ground state breaking the rotational symmetry and I φ,0 = 0 otherwise. Furthermore, he proved the inequality E (ω) φ,n > E (0) φ,n for any finite value of ω = 0. For the system in the thermal equilibrium at inverse temperature β, the same conclusion holds for the free energy, i.e., F (ω) φ,β > F (0) φ,β for ω = 0. These two inequalities consist of the "no-go theorem", which strictly prohibits the existence of spontaneously rotating time crystals for the ground state and thermal equilibrium state. Bruno also pointed out that the moment of inertia for an excited state I φ,n (n > 0) may be negative, which allows for the existence of time crystal at some excited eigenstate [11,12] or time-crystal Sisyphus dynamics near the energy minimum [13].
However, as Bruno's proof is restricted to a special model Hamiltonian, Wilczek and others later proposed new types of time crystals to avoid Bruno's no-go theorem [14,15]. A more general no-go theorem was given by Watanabe and Oshikawa (WO) [16]. Using the local order parameterφ( x, t), WO presented a precise definition for the time crystals, that is, the system is a time crystal if the correlation function lim V →∞ e iĤtφ ( x, 0)e −iĤtφ (0, 0) → f (t) (3) for large | x| (much greater than any microscopic scales) exhibits a nontrivial periodic oscillation in time t. In terms of the integrated order parameterΦ ≡ V d d xφ( x, 0), the condition reads equivalently lim V →∞ e iĤtΦ e −iĤtΦ /V 2 → f (t). At zero temperature, WO was able to prove the following equality for the ground state 1 where E 0 is the ground state energy and C is a constant dependent onΦ and H but not on t and V . The above inequality holds for the interaction which decays exponentially as a function of the distance or decays as power-law r −α (α > 0). Therefore, for any fixed time t, the correlation function f (t) remains constant in the limit V → ∞. The conclusion can be extended to the case of canonical ensemble at finite temperature by the help of Lieb-Robinson bound [17]. A related definition of time translation symmetry breaking (TTSB) was discussed in the language of C * algebras [18], which recovers WO's result for the absence of TTSB in equilibrium systems. It is still an open question that if one can bypass the limitations imposed by WO's no-go theorem and reach the desirable phase of a quantum time crystal in the ground state of an isolated quantum system without a driver [19].

Discrete time-translation symmetry breaking
The no-go theorem valid for equilibrium state does not exclude the possibility of spontaneous TTSB in nonequilibrium states, e.g., the intrinsically nonequilibrium setting of a periodically driven (Floquet) system, whose Hamiltonian has the discrete time-translation symmetry (DTTS) such that H(t) = H(t + T d ) with T d = 2π/ω d the period of driving field. It was found theoretically that the DTTS can be spontaneously broken in periodically driven atomic systems [20] and spin chains [21][22][23][24]. The DTTS is said to be broken if for every "short-range correlated" state |ψ , there exists an operatorΦ and integer n > 1 such that namely, the system responds at a fraction ω d /n of the original driving frequency [22].
For the detailed description of Floquet time crystals, it needs some basics for the periodically driven systems, i.e., the Floquet theory. In the next section, we will first introduce the Floquet theory and Magnus expansion, then describe in detail how to realise the Floquet time crystals by breaking DTTS.

Floquet theory
The DTTS of Floquet systems with periodic Hamiltonian H(t) = H(t + T ) enables us to use the Floquet formalism [34][35][36][37][38]. We write the Schrödinger equation for the quantum system as ( = 1) where x represents coordinates or spins. Similar to the Bloch theorem in solid state physics, the time periodicity of operator H (x, t) = H (x, t + T ) enforces the eigensolutions of Eq. (7) in the form of (so-called Floquet solution [34][35][36][37][38]) Here, |ϕ α (x, t) is called the Floquet mode and α is called the quasienergy, which is a real-valued variable modulo Ω = 2π/T . The terminology quasienergy is analogous to the quasimomentum characterizing the Bloch eigenstates in a crystalline solid. By substituting Eq. (8) into Eq. (7), we find According to Eq. (8), the new state |ϕ α,n (x, t) ≡ e inΩt |ϕ α (x, t) (n = 0, ±1, ±2, · · ·) is also the Floquet mode but with the shifted quasienergy α,n ≡ α + nΩ. Therefore, we can map all the quasienergies into a first Brillouin zone, e.g., 0 ≤ α < Ω. In order to calculate the quasienergy spectrum and the Floquet eigenstates of H (x, t), it is convenient to introduce a composite Hilbert space R ⊗ T , where R is the spatial (or spin) space and T is the space of functions with time periodicity T . We can choose the time-independent states |φ m (x) satisfying φ n (x)|φ m (x) = δ nm as the orthonormal basis of space R and the time-dependent Fourier vectors |e iqΩt with q = 0, ±1, ±2, . . . , as the orthonormal basis of space T . Then, we can obtain the matrix of H (x, t) in the complete set of {|φ m (x) ⊗ |e iqΩt } and calculate its quasienergy spectrum and the Floquet eigenstates.
On the one side, according to Eq. (8), the dynamics of the Floquet eigenstate |ψ α (x, t) at stroboscopic time moment t = nT (n ∈ Z) is given by On the other side, from the Schrödinger equation, the time evolution of the Floquet eigenstate |ψ α (x, t) can be directly written in the Dyson form where T is the time-ordering operator. Comparing Eq. (10) to Eq. (11), the stroboscopic dynamics of Floquet system can be described by an effective timeindependent Hamiltonian H F determined from or written formally as which has the eigenvalues α and eigenstates |ϕ α (x, 0) . The time-independent Hamiltonian H F is called Floquet Hamiltonian.

Magnus expansion
It is tempting to calculate the explicit form of the Floquet Hamiltonian (13). However, except very few examples like a driven harmonic oscillator, it is impossible to obtain a closed form of H F in general. In the regime where the driving frequency is much larger than the characteristic frequency of the system, the Floquet Hamiltonian can be given by the so-called Floquet-Magnus expansion. Given the time differential equation of motion for the time evolution operator i∂ t U (t, t 0 ) = H(t)U (t, t 0 ), the Magnus theorem [39,40] claims that the time evolution operator can be written in the form of U (t, t 0 ) = exp[−iΩ(t, t 0 )] with the operator Ω(t, t 0 ) determined by the following differential equation Bn n! x n . The first few nonzero Bernoulli numbers are B 0 = 1, B 1 = −1/2, B 2 = 1/6, B 4 = −1/30 and B 2m+1 = 0 for m ≥ 1.
We now replace H by εH in Eq. (14) and try a solution in the form of the Magnus series By substituting the above series (15) into Eq. (14) and equating powers of ε, we can obtain in principle all the expansion orders Ω n (t). For example, the first three orders are given by [40] , , For the explicit expressions of higher orders, one can find in Refs. [41][42][43][44][45]. We apply the above Magnus theorem to Floquet system with periodic Hamiltonian H(t) in the form of Fourier harmonics with the first three orders given by [46] H (0) The leading term H We should emphasise that, although the Magnus expansion method has been used for a long history in physics, the precise radius of convergence of Magnus series is still an open mathematical problem. A sufficient condition for the convergence of Magnus series (15) in the time interval t ∈ [0, t c ] is given by [40] tc 0 where the norm A is defined as the square root of the largest eigenvalue of A † A, and the operator H(t) needs to be a bounded operator in Hilbert space. More precise criterion of convergence can be found in Ref. [47]. For the Flouqet-Magnus expansion, the convergence condition becomes T 0 H(t) dt < π with T = 2π/Ω the time period. However, this convergence condition is not guaranteed in physical systems. For instance, since the generic nonintegrable infinite-size systems eventually heat up to infinite temperature for an arbitrary driving frequency, it was concluded that the radius of convergence for the Floquet-Magnus expansion vanishes in the thermodynamic limit and the quantum ergodicity results in the divergence of the expansion [48][49][50][51][52][53]. It is still unclear whether the divergence of the Floquet-Magnus expansion necessarily results in the heat up to infinite temperature due to the difficulty to estimate the radius of convergence for nonintegrable many-body systems. It was recently conjectured that the divergence of the Floquet-Magnus expansion is a universal feature of periodically driven nonlinear systems with energy localization [54]. Nevertheless, since the system costs exponentially long time to reach the infinitetemperature featureless state if the driving frequency is much larger than the local energy scale, the prethermal dynamics can be well described by a truncation of the Floquet-Magnus expansion [53, 55-61].

Subharmonic modes
To understand the formation of Floquet time crystals, let us consider a one dimensional chain of N spin-half particles with a binary driving protocol where σ i are Pauli operators, B z i ∈ [0, W ] is a random longitudinal field and δ is the imperfection of the uniform transversal field B x . Here, we set the dimensionless time interval τ = π n with n ∈ N and T d > τ /B x is the driving period. For the case of perfect driving (δ = 0) and clean system (disorder strength W = 0), the Floquet unitary (12) has the form ( = 1) where we have chosen the Floquet period T = T d . Obviously, the short-range correlated state |φ 0 ≡ | ↑↑↑ · · · ↑ is NOT the eigenstate of U f since the operator R τ rotates all the spins. The real Floquet eigenstates |ψ m (m = 0, 1, · · · , n − 1) can be constructed in the following form where |φ q ≡ R q τ | ↑↑↑ · · · ↑ (q = 0, 1, · · · , n − 1). The Floquet eigenstates (21) have Z n symmetry, i.e., R n τ |ψ m = |ψ m from R τ |ψ m = e −i2mτ |ψ m . Since the interaction part in H 2 is invariant by the unitary transformation of operator R τ , we have the quasienergy of Floquet eigenstate |ψ m from U f |ψ m = e −i m T d |ψ m where ω d = 2π/T d is the driving frequency. The set of n Floquet eigenstates |ψ m (m = 0, 1, · · · , n − 1) form a multiplet in the quasienergy spectrum with equidistant quantized quasienergy ω d /n as illustrated in Fig. 1(left). Note that the short-ranged correlated basis |φ q in Eq. (21) are not orthogonal for a finite-size system, i.e., φ q |φ q = cos (q−q )π n N . However, in the thermodynamic limit N → ∞, all the short-ranged correlated states |φ q are orthogonal φ q |φ q = lim N →∞ cos (q−q )π n N → δ qq . The long-range correlated eigenstates |ψ m are more like multi-component Schrödinger-cat states, which are unlike to survive in the thermodynamic limit. If the system is initially prepared in one Floquet eigenstate |ψ m , any local perturbation in the environment can induce the system to collapse randomly into one of the short-range correlated basis |φ q . As consequence, the Z n symmetry of |ψ m is spontaneously broken in the thermodynamic limit. From Eq. (21), we can express the short-range correlated states with Floquet eigenstates e −i2mqτ |ψ m , with q = 0, 1, · · · , n − 1.
Using Eq. (22), we have the stroboscopic time evolution of |φ q as following We see that φ q (nT d )|Φ|φ q (nT d ) = φ q (0)|Φ|φ q (0) with integer n > 1. According to the definition Eq. (5), the DTTS is spontaneously broken in the thermodynamic limit accompanied by the breaking of a global internal Z n symmetry described by R τ . The symmetry-broken localised states are also called subharmonic modes since their oscillating frequency is a fraction ω d /n of the driving frequency ω d . The Z n symmetry breaking (period-n tupling) was discussed, e.g., in the n-hands clock model [62,63], ultracold bosons in a ring lattice [64] and long-range interacting systems [65]. Since the original time-dependent Hamiltonian also satisfies H(t) = H(t + nT d ), we have the freedom to choose the Floquet period T = nT d . In this setting, the Floquet Hamiltonian is given formally by Under this choice, the n Floquet multiplets are "degenerate" within the Floquet-Brillouin zone 0 ≤ α < ω d /n as illustrated by the black lines in Fig. 1(middle). They are long-range correlated eigenstates analogous to multi-component Schrödinger-cat states with Z n symmetry as shown by Eq. (21). In the thermodynamic limit, the n long-range correlated states randomly collapse into one of the n short-range correlated states, which are illustrated by the coloured wave packets in Fig. 1(middle). As consequence, the Z n symmetry of long-range correlated eigenstates is spontaneously broken, and the resultant n short-range correlated states exhibit oscillations with the same time period nT but with phases differed by 2π/n [see Eqs. (23) and (24)]. We arrange the n short-range correlated states on a Z n circle as shown in Fig. 1(middle). The emergence of a multiplet of Floquet eigenstates with equidistant quasienergy ω d /n as illustrated in Fig. 1(left) is the prerequisite for DTTS breaking [21][22][23].

Stabilising mechanisms
However, the subharmonic modes do not necessarily mean the Floquet time crystals.
In the spin chain model discussed above, the subharmonic modes are generally unstable since the system is generically heated up to a featureless infinite-temperature state by the driving [49][50][51]. As a solution, disorder [e.g., setting W > 0 for random field B z i in Eq. (19)] is used to stabilise the subharmonic modes [21][22][23][24] via the mechanism of many-body localization (MBL) [66][67][68][69]. To determine if the system is in the MBL phase, one can calculate the quasienergy statistics ratio, r = min(∆ m , ∆ m+1 )/ max(∆ m , ∆ m+1 ), where ∆ m = m+1 − m is the mth quasienergy gap [21,70], by averaging over both disorder and the quasienergy spectrum. In the thermal phase, this value approaches the circular orthogonal ensemble of r ≈ 0.527 while, in the MBL phase, the value should approach the Poisson limit of r ≈ 0.386 [50]. If the system is in the MBL phase, it was shown that [24] the spontaneous breaking of DTTS is rigid for a certain range of driving imperfection 0 < δ < δ c . If the imperfection is larger than the critical point δ > δ c , the system is still in the MBL phase but the DTTS is not broken. Floquet time crystals can also exist in clean Floquet systems without disorder [27,28,65,[71][72][73][74]. Although the fate of a generic driven many-body system is a trivial infinite temperature state with no long-range order, there exists a prethermal state with an extremely long lifetime provided that the driving frequency is much larger than the local energy scales in the system [56,57,59,75]. The system needs many local rearrangements of the degrees of freedom to absorb an "energy quanta" from the driving, which is a high-order process taking exponentially long time even in a clean system with no MBL and without disorder. The DTTS broken state in this long time prethermal process is called prethermal time crystal [76,77].
By coupling the Floquet many-body system to a cold bath, it is also possible to protect the prethermal time crystals for infinitely-long time [76]. Actually, the spontaneous breaking of DTTS was already discussed even before Wilczek's seminal work in a modulated dissipative cold atom system [78], and the DTTS symmetry breaking has been reported in the experiment [79]. However, strictly speaking, Floquet time crystals by our definition are the quantum phenomena emerging in closed systems not in open systems. Therefore, we exclude the "time order" or the self-organised phenomena in driven-dissipative systems introduced by Prigogine [80], such as synchronisation [81] and Belousov-Zhabotinsky reaction [82], from the family of Floquet time crystals.
In fact, MBL with disorder spin systems is NOT a necessary condition to realise Floquet time crystals in other physical systems (e.g., the quantum gas). There exist other mechanisms to break Z n symmetry and stabilise the broken states other than MBL. Imagine the Floquet Hamiltonian (25) has a Z n -lattice potential along the Z n circle as illustrated by the blue curve in Fig. 1(right). If the lattice potential is strong enough, the Floquet eigenstates are the superposition of localised states near the lattice sites (Wannier functions), which can tunnel to their neighbouring sites with hopping rate J. In this case, the system can be described by a tight-binding model. As a result, the quasienergies of the Floquet eigenstates form a band (∝ |J|) as indicated in Fig. 1(right). The tunnelling rate can be tuned by system parameters. In the limit J → 0, the Floquet eigenstates become degenerate and form a multiplet ready for symmetry breaking. At this point, any nonequilibrium fluctuations can induce the Floquet eigenstates to collapse into one of the localised states, and the broken state is protected by the local potential well of Z n lattice. If the subharmonic modes have attractive interaction, the spontaneous symmetry breaking occurs when the interaction strength is larger than a critical value proportional to the tunnelling rate J [20]. This mechanism of realizing time crystals without referring to MBL has also been studied in the driven Lipkin-Meshkov-Glick (LMG) model [71], where only the collective degree of freedom of spins is relevant due to conservation laws and in the thermodynamic limit the collective spins behave classically, therefore subharmonic response is protected by nonlinear resonances. Time crystals can also be protected by Floquet dynamical symmetry (FDS) and robust against perturbations respecting the FDS symmetry [83].

Condensed matter physics in time crystals
Most works so far on time crystals focus on the spontaneous breaking process of DTTS and the possible mechanisms to stabilise the subharmonic mode. However, the interplay of subharmonic modes is still not fully explored. We call the manybody physics of subharmonic modes the field of condensed matter physics in time crystals. In this section, we review the very preliminary results in this field based on the models of periodically driven quantum oscillators (e.g., driven quantum gas), instead of periodically driven spin models since there is no such work in spin models yet. In general, the Hamiltonian of driven quantum oscillators in one dimension is described by where H s (x i ,p i , t) is the single-particle Hamiltonian which is NOT assumed to be periodic in space but periodically time depends on the distance of two particles. We first discuss the physics of single-particle Hamiltonian from Sec. 3.1 to Sec. 3.4, and then include the effects of interaction from Sec. 3.5 to Sec. 3.8.

Phase space lattices
When talking about condensed matter physics, there is usually a lattice symmetry assumed in the Hamiltonian. However, Hamiltonian (26) has no such lattice structure in real space, and the time periodicity is also not the lattice structure we refer to since this time periodicity is already included in the Floquet theory. In fact, the lattice structure for the condensed matter physics in Floquet time crystals is hidden in phase space. The general idea goes as follows [84]. Let us begin with a classical singleparticle Hamiltonian described by a time-independent H 0 (x, p). Since the energy is conserved in this case, the classical orbit is just the isoenergetic contour of H 0 (x, p).
We define the action variable by along the isoenergetic contour of H 0 (x, p) and rewrite the Hamiltonian as a function of action H 0 (I). The conjugate variable of action I is called angle θ, which does not explicitly appear in H 0 (I). According to the canonical equation of motion, we have θ(t) = Ωt + θ(0) with Ω(I) = dH 0 (I)/dI representing the frequency of the regular motion orbit. We further introduce a perturbation term of the form and the spatial function can be expressed in terms of the action and angle variables h(x) = l h l (I)e ilθ . Thus the whole single-particle Hamiltonian is given by, Now if we tune the driving frequency to meet the resonant condition ω d = nΩ(I 0 ) with n ∈ N and I 0 the action of a resonant orbit, in the vicinity of this resonant trajectory the motion of the particle can be separated into a fast oscillating and a slow varying parts, where the latter can be described by a new angle variable Θ = θ − ω d t/n. Using the generating function G 2 = I(θ − ωt/n), we transform into the rotating frame with frequency ω d /n and retain the time-independent Hamiltonian from Eq. (28) in the rotating wave approximation (RWA) where V eff (I, Θ) = k h kn (I)f −k e iknΘ is a periodic potential in the angel direction Θ of rotating frame. This effective Hamiltonian essentially describes a particle moving on the Θ ring in the presence of a time-independent effective lattice potential V eff (I, Θ) which can be engineered by tuning the Fourier components f k and h l in the experiments. For example, for the monochromatic driving with frequency ω d , the effective lattice potential is V eff (I, Θ) ∝ cos(nΘ) meaning the driving resonates with the classical orbit with frequency ω d /n and generates n identical local wells uniformly distributed on the Θ ring. Below, we will introduce several concrete examples to show various lattice structures in phase space. The first example is a nonlinear oscillator driven by an external field coupling nonlinearly to the coordinate, with the single-particle Hamiltonian reading [85] Here ω 0 is the frequency of the oscillator and ω d the driving frequency. The nonlinearity of coupling is characterized by the exponent n. If n = 1 the model T / Td (30) is the linearly driven Duffing oscillator [86][87][88][89], for n = 2 it is a parametrically driven oscillator [90][91][92]. The model (30) extends the exponent n to be any integer larger than two, and has been recently realised up to n = 5 in the experiment with superconducting circuits [93]. We are interested in the subresonance regime, where the driving frequency ω d is close to n times ω 0 , i.e., the detuning δω ≡ ω 0 − ω d /n is much smaller than ω 0 . By taking the Harmonic oscillator as H 0 (x, p), the time-independent RWA Hamiltonian with action and angle variables is given by Here, the Hamiltonian has been scaled by energy unit 4m 2 ω 2 0 δω 2 /3ν and µ is the scaled driving strength given by µ = 3νf 2 mω0δω −3ν with red detuning δω < 0.
The RWA Hamiltonian (31) displays a new symmetry not visible in (30). To show the lattice structure, we define two alternative conjugate variables and plot H RW A in phase space spanned by X and P in Fig. 2(a). Clearly, there are n identical local wells distributed uniformly on a circle with radius √ 2I 0 with I 0 = 1/2. The Hamiltonian shows a Z n lattice structure in phase space. The n local wells only exist in a certain range of the driving strength 0 < |µ| < µ c , where the critical driving strength is given by [85]. The second example is a harmonic oscillator in driven optical lattice, where the power-law driving in model (30) is replaced by a cosine-type driving, i.e., the single- Square lattice Hexagonal lattice Quasicrystal particle Hamiltonian is given by [94] H s (t) = The above Hamiltonian also describes a cavity mode driven by dc-biased Josephson junction(s) [95][96][97], which has triggered a new research field of Josephson photonics [98][99][100][101]. Under the subresonance condition ω d = nω, the RWA Hamiltonian in the classical limit with action and angle variables is given by [94] H RW A (I, Θ) = 2µJ n ( where J n (•) is the Bessel function of order n. The Hamiltonian has been scaled by energy unit mω 2 /k 2 and µ ≡ Ak 2 /mω 2 is the dimensionless driving strength. In Fig. 2(b), we plot the scaled RWA Hamiltonian in phase space. The whole lattice structure clearly shows a Z n symmetry. The zero lines (i.e., H RW A = 0) divide the whole phase space lattice into many small " cells ". The center of each cell is a stable point corresponding to either a local minimum or a local maximum. The area inside the cell represents the basin of attraction for the stable state in the center. The RWA Hamiltonian shows a multiple Z n lattice structure in phase space.
In the above two examples, the periodic structure is actually one dimensional, i.e., the lattice structure appears along the angular Θ direction in phase space. One may ask if it is possible to create a two dimensional lattice in the whole phase space. The answer is yes [102]. Let us consider the famous model of kicked harmonic oscillators (KHO) [103][104][105][106] with the single-particle Hamiltonian given by where K is the dimensionless kicking strength and τ is the dimensionless kicking period. In the resonant condition that the kicking period satisfies τ = 2π/q 0 with q 0 an integer and the weak kicking strength limit|K| 1, the single-particle dynamics is still dominated by the fast harmonic oscillation but with slowly modulated amplitude and phase. The RWA part of the single-particle Hamiltonian is given by [102] H RWA (X, P ) = K q 0 q0 j=1 cos X cos 2πj In Fig. 3, we plot H RWA (X, P ) in phase space for different q 0 . We see that the H RWA (X, P ) has a square lattice structure for q 0 = 4, hexagonal lattice structure for q 0 = 3 or q 0 = 6, and even quasicrystal structure for q 0 = 5 or q 0 ≥ 7. There is another way to view the hidden lattice structure of Hamiltonian. In the laboratory frame, the crystalline structure in Θ is reproduced in the time domain as the relation between Θ and t is linear, i.e., Θ = θ − ω d t/n. Imagine we measure the probability of the system by locating a detector at a fixed spatial point close to the classical resonant trajectory, e.g., at point (1, 0) in Fig. 2(a). If the system is in the subharmonic mode, i.e., a localised wave packet at one site of the Z n lattice, the detector will be triggered with longer time period nT d breaking the original DTTS of driving field as indicated by the lower subfigure in Fig. 2(c). If the system is in one Floquet eigenstate, which is the superposition of n localised states in the Z n lattice, the detector will respond with the driving period T d as indicated by the upper subfigure in Fig. 2(c). This time-domain lattice picture was introduced by Sacha [107] to study cold atoms bouncing on an oscillating mirror. In this picture, the roles of space and time are exchanged making it possible to discuss analogous condensed matter phenomena like Anderson localisation and Mott insulator in the time domain. Here in this review, we keep our discussion in the rotating frame and adopt the picture of phase space lattices. In quantum dynamics, one has to find a unitary operator to describe the discrete lattice transformation under which the RWA Hamiltonian is invariant. Such operator is natural to be defined in the picture of phase space lattice, and allows us to develop a band theory, which is to be discussed in the next section.

Band theory
In the above section, we discuss the phase space lattices in the classical limit. In quantum mechanics, the lattice symmetry results in the band structure of the energy spectrum. Let us start from the example Hamiltonian (30). We transform into the rotating frame via the unitary operatorÛ = e i(ω d /n)â †â t , whereâ is the annihilation operator of the oscillator, resulting in The RWA part Hamiltonian is given by (scaled by energy unit of 4m 2 ω 2 0 δω 2 /3ν) Here, λ ≡ −3ν /(4m 2 ω 2 0 δω) is the scaled dimensionless Planck constant, which describes the quantumness of our system. In the classical limit λ → 0, expression (37) goes back to the classical form (31). The non-RWA part has the following explicit form [85] H non−RW A = 1 3 According to the Floquet-Magnus expansion (18), the RWA HamiltonianĤ RW A is just the leading order of Floquet Hamiltonian, i.e.,Ĥ RW A = H (0) The RWA Hamiltonian (37) displays a new symmetry not visible in the original Hamiltonian (30). We define a unitary operatorT τ ≡ e −iτâ †â , which has the propertieŝ T † τ aT τ = ae −iτ andT † τ a nT τ = a n e −inτ . It is direct to find that the RWA Hamiltonian is invariant under discrete transformation According to the Bloch's theorem, the eigenstates ψ m (Θ) of the RWA Hamiltonian, H RW A ψ m (Θ) = g(m)ψ m (Θ), have the form ψ m (Θ) = ϕ m (Θ)e −imΘ , with a periodic function ϕ m (Θ + τ ) = ϕ m (Θ). Here, the integer number m, which is called quasinumber in Ref. [85], plays the role of the quasi-momentum k in a crystal. While the quasi-momentum k is conjugate to the coordinate, the quasi-number m is conjugate to the phase Θ. In Fig. 4(a), we plot the quasienergy band structure in the reduced Brillouin zone mτ ∈ (−π, π] from exact numerical diagonalization. Here we relabel the eigenstates ψ m (Θ) by ψ lm (Θ), where l = 1, 2, ... is the label of the bands counted from the bottom. For finite values of n (in our numerical simulation we chose n = 10) the quasienergy band spectrum is discrete. It would become more continuous in the limit of large n. For the multiple Z n lattice Hamiltonian (33), a similar quasi-number band theory was developed in Ref. [94].

Tight-binding model
The phase space lattice is fundamentally different from the real space lattice as the two degrees of freedom of phase space do not commute, i.e., [X,P ] = iλ for the Hamiltonian (37). As a result, the band structure of a phase space lattice shows some novel features. For the band structure shown in Fig. 4(a), the spectrum of the lower l-th band is approximately given by a tight-binding model [85] g l (m) = E l − 2|J l | cos(mτ + δτ ), centered around E l and with bandwidth d l = 4|J l |. The result shows a surprising asymmetry. From the plot of the phase space lattice in Fig. 2(a) ) with two complex numbers ξ 1 , ξ 2 determining the moving path. The geometric phase factor is given by e i 1 λ S with S = 1 2 Im[ξ 2 ξ * 1 ] being the area of the enclosed path (red area). expected a degeneracy g(m) = g(−m), since clockwise and anticlockwise motion should be equivalent, as in the case of orbital motion. However, in the noncommutative phase space, the concept of point is meaningless. Instead, we should define a coherent state |α which is the eigenstate of the lowering operator, i.e.,â|α = α|α witĥ a ≡ (X + iP )/ √ 2λ. As shown in Fig. 5, we observe that a coherent state moving along a closed path in phase space naturally acquires an additional quantum phase factor e iS/λ , where S is the enclosed area [108]. This geometric phase is responsible for the asymmetry of quasienergy band structure, and also called Peierls phase [109] for describing the hopping of charged particles between tight-binding lattice sites in a magnetic field. For neutral ultracold atoms in an optical lattice, a controlled Peierls phase can be created by shaking the lattice to realize artificial gauge fields [110][111][112].
According to the Bloch theorem, the tight-binding eigenstate is where |φ q,l is the l-th localised Wannier level in the q-th well of the Z n lattice, and T τ = e −iτâ †â is the discrete rotation operator shown in Eq. (39). Reminiscent of eigenstate expression (21) for the spin model, they both reflect the underlying Z n symmetry of corresponding systems. The nearest coupling in Eq. (40) is given by where we have written the Wannier functions in the angle Θ representation. Interestingly, the calculation in Ref. [85] shows the nearest neighbor coupling rate is in general a complex number J l = |J l |e iτr 2 /2λ wherer is the average radius. Hence the asymmetry factor is δ =r 2 /2λ (mod n). In Fig. 4(c), we show the dependence of the asymmetry factor δ on the driving strength µ, obtained both in the tight-binding calculation described above and from a numerical simulation. The asymmetry factor δ =r 2 /2λ (mod n) in Eq. (40) from the imaginary part of J l is valid for large exponent n 1. More precise calculation for small exponent (period tripling n = 3) was given in Ref. [113]. The complex tunnelling parameter J l naturally arises in the plane of the phase space due to the noncommutative geometry.
The amplitude of tunnelling rate J l was calculated using the WKB semiclassical approximation in Ref. [85], i.e., |J l | = λωe 2π exp − rs λ W l with the exponent W l depending on the system parameters. Here, r s is the radius of the saddle points of H RW A and ω e is the harmonic frequency near the local minimal points. In Fig. 4(b) we compare our approximate result for the first bandwidth d 1 = 4|J 1 | versus driving µ to numerical results. In the tight binding regime they agree well with each other. The tunnelling rate |J l | decreases exponentially to zero in the classical limit of λ → 0. This means the n eigenstates in each quasienergy band become degenerate and form the multiplet in the Floquet spectrum as shown in Fig. 1(middle). In this limit, the Floquet eigenstate with Z n symmetry cannot survive and will be broken into one of the localised Wannier states. Therefore, Z n symmetry is spontaneously broken and the broken state oscillates with time period of nT , i.e., the subharmonic mode of Floquet time crystal. This period-multiplication phenomena predicted by the model (30) has been verified in the experiment with superconducting circuits [93,114], where experimentalists observed robust output radiation (subharmonic oscillations [115]) at a frequency close to a fraction ω d /n of the driving frequency and evenly 2π/n-shifted phase components. These observations call for further exploration of the quantum properties of the period-multiplication and a search for engineering multicomponent quantum superpositions of macroscopic coherent states (multicomponent cat-states). While most of these experimental findings display classical fixed points and the hopping between them induced by classical noise, the quantum mechanical coherences in the form of multi-component cat-states are still absent. The Z 2 symmetry breaking of parametrically driven model has been reported in the experiment with cold atoms [79] and analysed in Ref. [78]. Recently, there is a lot of interest in Z 3 -lattice of triply driven oscillator [113,114,[116][117][118], e.g., nonlocal random walk [116], quantum state preparation [117] and dissipative phase transition [118]. A similar calculation for the band structure of the multiple Z n lattice (34) with tight-binding model has been provided in Ref. [94]. There, the situation is more complicated since the tunnelling between nearest Z n lattices also play an important role when bands cross each other.
In the above discussion, we work in the single-particle picture. For an ensemble of identical particles free of interaction, we can describe the corresponding many-body Hamiltonian with the following tight-bing model Here,â q,l is the annihilation operator of the Wannier state |φ q,l . We have restricted to the Fock space |n 0,l , n 1,l , · · · , n q=n−1,l with n q,l the particle number occupying the Wannier state |φ q,l . A similar tight-binding model was introduced independently by K. Sacha in the study of bouncing atoms on an oscillating mirror with the subresonant codition [107].

Topology
Topological phenomena in condensed matter physics have been intensively investigated in the past decades. In quantum Hall physics, the origin of topology comes from the geometric phase induced by the applied magnetic field and the resultant band can be characterized by its topological invariant (Chern number or TKNN invariant), which relates to the quantized Hall conductance directly. In the noncommutative phase space, a geometric phase factor of a quantum state appears naturally when it moves along an enclosed path as shown in Fig. 5. We would expect that similar topological phenomena also exist in phase space lattices. This is actually true for the example Hamiltonian (36). In the case of square lattice (q 0 = 4), Hamiltonian (36) is further simplified as This Hamiltonian is closely related to the established Harper's equation, which is a tight-binding model governing the motion of noninteracting electrons in the presence of a two-dimensional periodic potential and a uniformly threading magnetic field [109,119]. The H sq (X,P ) is invariant under discrete translation in phase space by two operatorsT 1 ≡ e i 2π λP andT 2 ≡ e i 2π λX . However, the translation operatorsT 1 and T 2 are NOT commutative in general, i.e., with integer powers r, s ∈ Z except for 2πrs/λ ∈ Z. If λ/2π = p/q with p and q two coprime integers, we choose the following two commutative translational operators (r = 1, s = p):T X ≡T 1 = e i 2π λP ,T P ≡T p 2 = e i 2πp λX . Therefore, we can search for the common eigenstates of the commutative operatorsT X andT P with eigenvalues given by e i2πk X and e i2πpk P respectively. The boundaries of the two dimensional Brillouin zone are defined by 0 ≤ k X ≤ 1 and 0 ≤ k P ≤ 1/p, where k X and k P are quasimomentum and quasicoordinate, respectively [120]. In Fig. 6(a), we plot the quasienergy spectrum for λ/2π ∈ [0, 1], showing a Hofstadter's butterfly structure identical to that in quantum Hall systems. In Fig. 6(b), (d) and (e), we plot the quasienergy band structures in the two-dimensional Brillouin zone (k X , k p ) for λ/2π = 1/2, 1/3 and 2/3, respectively. We see that the quasienergy band structure is two-fold degenerate for λ/2π = 2/3 while there is no degeneracy for λ/2π = 1/2 and λ/2π = 1/3. In fact, for each rational λ/2π = p/q (remembering p, q are coprime integers), the spectrum contains q bands and each band has a p-fold degeneracy [121]. To show the underlying topology of the quasienergy bands, one can calculate the Chern number for a given band [122] where the contour C is integrated over the boundary of the Brillouin zone. Accordingly, the Chern number of a gap below (above) the zero energy line can be defined as the sum of Chern numbers of all the quasienergy bands below (above) the gap [102]. The Chern number of some gaps are calculated and labelled in Fig. 6(a).
The topology manifests itself also in the dynamics. In reality, the driven oscillators are inevitably in contact with the environment resulting in dissipation or decoherence of the quantum system. In the Born-Markovian approximation, the dissipative dynamics of the quantum KHO is described by the following master equation [102], where the Lindblad superoperator is defined by D[Ô]ρ ≡ÔρÔ † − 1 2 (Ô †Ô ρ + ρÔ †Ô ) and κ characterizes the dissipation rate and n 0 is the Bose-Einstein distribution of the thermal bath. It should be emphasized that the above dissipative dynamics of KHO is based on the original full Hamiltonian (35) without RWA. To visualize the time evolution of quantum state, we define the Husimi Q-function of a given density matrix [123] where |α is the coherent state. In Ref. [102], the time evolution of the Husimi Qfunction starting from the ground state of the harmonic oscillator was calculated, which is shown in Fig. 7. It is clearly seen that a final steady state with square lattice structure in phase space forms gradually revealing the underlying square structure of Hamiltonian (44). Interestingly, the transient state shown in Fig. 7(b) has no reflection symmetries with respect toX andP albeit the RWA Hamiltonian (44) has, which results in a chiral feature as marked by the dashed lines along the backbone of the quasiprobability distribution. This chirality is a reflection of the topological property of our system and the noncommutative geometry of the phase space. The probability amplitude of a particle appearing at a fixed point (a coherent state in phase space) is the sum of all the possible trajectories. Different from the path integral in 2D real space, each trajectory in phase space associates with a geometric phase due to the noncommutative geometry. The interference of geometric phases breaks the mirror symmetry of phase space. As approaching the stationary state, the chirality disappears in the end since the stationary state should recover the mirror symmetry of the RWA Hamiltonian (44). In Ref. [124], the Su-Schrieffer-Heeger (SSH) model [125] was proposed to be realised using bouncing atoms on an oscillating mirror. In the SSH model, the topological phase exhibits zero energy edge states protected by the topology of the bulk. In the same paper, the authors also proposed to simulate the Bose-Hubbard Hamiltonian with repulsive on-site and nearest-neighbor interactions, where the topological Haldane insulator phase [126] may exist. But in order to understand the strongly correlated phenomena, we need to provide a general way to deal with the interaction term in Hamiltonian (26), which is the centre topic in the next section.

Interactions
Until now, we have neglected the interaction terms in the original Hamiltonian (26). Including the interactions between particles, the total RWA Hamiltonian can be written as Here, H RW A (X i ,P i ) is the single-particle RWA Hamiltonian discussed above. The interaction term U (X i ,P i ;X j ,P j ) is the RWA part of the transformed real space interaction V (x i −x j ) in the rotating frame. In general, the interaction term U (X i ,P i ;X j ,P j ) is defined in the phase space of the rotating frame and depends on both coordinates and momenta of two particles. Below, we discuss the explicit form of U (X i ,P i ;X j ,P j ) both in the classical limit and in quantum mechanics.
In classical dynamics, we can omit the hats of operators in the Hamiltonian (46). For a given real space interaction potential V (x i − x j ), a general form of U (X i , P i ; X j , P j ) has been found in Ref. [127] , and J 0 (•) is the Bessel function of zeroth order. The second equality comes from the integral representation of the Bessel function, i.e., J 0 (x) = (2π) −1 π −π exp(−ix sin τ )dτ , and shows that U (R ij ) is in fact the time average of the original interaction V (x i − x j ) over the oscillation period. Since U (R ij ) is defined in the phase space of the rotating frame and it is only a function of the phase space distance R ij , we call it the phase space interaction. However, there is a divergence problem in calculating the phase space interaction U (R ij ) for, e.g., the Coulomb interaction. In Ref. [127], the origin of divergence was analysed and the correct U (R ij ) was obtained by introducing a renormalization procedure. For instance, the phase space interaction of the contact interaction in real space , which is used to describe the effective interaction between neutral ultra cold atoms, is given by U (R ij ) = π −1 β/R ij . A more general interaction potential form is the power-law interaction potential, i.e., V (x i − x j ) = β 2n /|x i − x j | 2n with integer and half integer n ≥ 1/2. If n = 1/2, the potential is the Coulomb potential. If n → ∞, the potential becomes the hardsphere potential with a diameter β. The renormalized phase space interaction for this power-law interaction is for n = 1, 3 2 , 2, 5 2 , · · · β π R ij , for n → ∞.
Here, γ is the collision factor given by γ = (4n − 1) 1 2n−1 . The real space power-law interaction is plotted in Fig. 8(a) and the corresponding phase space interaction is plotted in Fig. 8(b). Several interesting features of phase space interaction should be pointed out: (1) for the Coulomb potential (n = 1/2), U (R ij ) still keeps the form of Coulomb's law, up to a logarithmic prefactor (renormalised "charge" β); (2) for the special case of n = 1, U (R ij ) is a constant meaning there is no effective interaction in phase space; (3) for the special case of n > 1, U (R ij ) even grows with R ij ; (4) for the hard-sphere interaction (n → ∞), U (R ij ) increases linearly with phase space distance, reminiscent of the confinement interaction between quarks in QCD. All these predictions have been justified by simulating the many-body dynamics numerically in Refs. [102,127].
In quantum mechanics, the phase space interaction introduced above for the classical dynamics becomes subtle since the phase space is a noncommutative space.  As discussed in Fig. 5, the concept of point is meaningless in the noncommutative phase space. Instead, it is only allowed to define the coherent state |α as the point in noncommutative geometry. Similarly, the concept of distance also needs to be reexamined. It is natural to define the phase space distance operator bŷ From the noncommutative relationship [X i , where the operators ∆X ij ≡X i −X j , ∆P ij ≡P i −P j are the relative displacement of two particles in phase space. Reminiscent of the Hamiltonian operator of a harmonic oscillator, we see that the phase space distance is a quantized variable with a lower limit √ 2λ. In Ref. [102], a compact expression for the phase space interaction in the operator form was given by Here, V q = (2π) −1 +∞ −∞ dxV (x) exp(−iqx) is the Fourier coefficient of the real space interaction V (x i − x j ) and L n (•) is the Laguerre polynomials.
If the quantum states of two particles are two wave packet φ(i) and ϕ(j) as shown in Fig. 9(a), their averaged phase space interaction is given by The explicit expressions of U ij for the two displaced coherent states were calculated in Ref. [102]. For the contact interaction V (x 1 − x 2 ) = βδ(x 1 − x 2 ), the resultant U ij is a function of the phase space distance between the centers of two coherent states R ij , i.e., Here, I 0 (•) is the zeroth order modified Bessel function of the first kind. For the hard-sphere interaction, i.e., V (x i − x j ) = β 2n /|x i − x j | 2n for n → ∞, the result is In Fig. 9, the U (R ij ) as a function of R ij and its long-range asymptotic behaviour are plotted. We see that a point-like contact interaction indeed becomes a Coulomb-like interaction U (R ij ) ∼ π −1 β/R ij in the long-distance limit R ij 2 √ λ, and the phase space interaction of the hard-sphere potential becomes linear in the long-distance limit U (R ij ) → βπ −1 R ij , which are consistent with the pure classical analysis.

Mott insulator
With the theory of phase space interaction, we can extend the tight-binding model (43) for free particles to the many-body Hamiltonian with interactions. Assuming the particles are spinless Bosons, we can write the Floquet many-body Hamiltonian for the l-th band directly Here, J l is the tunnelling rate between nearest lattice sites denoted by i, j , and n i,l =â † i,lâ i,l is the number operator occupying the Fock space | · · · , n i,l , · · · in the Wannier representation. Model Hamiltonian (54) is valid if the interaction U ij,l is weak compared to the quasienergy gaps. If the phase space interaction U ij,l does not have compact form as given by Eq. (52) or it is difficult to obtain an analytical form, we can calculate it numerically from Eq. (51) or refer to the original expression in the laboratory frame where φ i,l (x m , t) represents the l-th Wannier wave function on the phase space lattice site i for the m-th particle in the laboratory frame. For the contact interaction V (x m − x n ) = g 0 δ(x m − x n ), we have U ij,l = g 0 T 0 dt dx|φ i,l (x, t)| 2 |φ j,l (x, t)| 2 . This is the method adopted by Sacha et al in their study [20,84,107,124,[128][129][130][131].
In Ref. [107], Sacha et al discussed the Mott transition using the model of bouncing atoms on an oscillating mirror. In the limit of free interaction, the ground state of many-body Hamiltonian (54) is a superfluid state with long-time phase coherence. However, for sufficiently strong repulsive contact interaction (g 0 > 0) such that U ii,l N |J l |/n (neglecting off-site interaction terms), the ground state of (54) is a single Fock state with well defined numbers of atoms occupying each Wannier state, i.e., |N/n, N/n, · · · , N/n with N denoting the total number of particles. As a result, a gap opens between the ground state level and excited levels, thus consequently a Mott insulator phase is achieved. Sacha et al also investigated the model (54) in two dimensional space by extending atoms bouncing between two orthogonal harmonically oscillating mirrors [130,132]. By including off-site interactions in (54), more intriguing condensed matter phenomena, e.g., the topological Haldane phase mentioned in the end of Sec. 3.4 , are possible to be explored [124].

Heisenberg model
If the two particles are indistinguishable and have spins (i.e., Fermions with half spinŝ i ), their spatial state is either antisymmetric or symmetric depending on the symmetry of total spin state, i.e., For singlet (triplet) spin state the wave function is symmetric (antisymmetric). Therefore, the average phase space interaction potential of the (anti)symmetric state is Here, the direct interaction part U ij has been given by Eq. (51) while the the exchange interaction part is The direct interaction part U ij corresponds to the classical interaction while the exchange interaction part J ij is a pure quantum effect without classical counterpart, which we call the Floquet exchange interaction energy. The Floquet exchange interaction energy of two coherent states was calculated from Eq. (55) in Ref. [102]. For the contact interaction, it is found that the exchange interaction is equal to the direct interaction, i.e., J(R ij ) = U (R ij ). The equality comes from the δ function modelling the contact interaction and the fact that the spatial antisymmetric state of two particles has zero probability to touch each other. If the interaction potential differs from the δ-function model, the phase space interaction U ij and the Floquet exchange interaction J ij can behave independently. For example, if the two particles are charged and confined in quasi-1D ion trap, the direct interaction and the Floquet exchange interaction were calculated in Ref. [133], i.e., U (R ij ) ∝ 1/R ij and J(R ij ) ∝ /R 3 ij . In this case, the exchange interaction decays in the form of inverse-cube law and is a quantum effect due to the prefactor . One should always remember that we are studying the stroboscopic dynamics of a periodically driven system. The two particles collide with each other during every stroboscopic time step, and the exchange interaction comes from the overlap of their wave functions during the collision process. For the hard-sphere interaction, the particles cannot cross each other [102], and thus the Floquet exchange interaction is zero J ij = 0.
If the particles are tightly bound at their equilibrium points by the phase space lattice, the direct phase space interaction U (R ij ) does not play a role in the dynamics. The dynamics of the spins is determined by the Floquet exchange interaction and is described by the Heisenberg model with isotropic spin-spin interaction, i.e., In Fig. 10, we show the Heisenberg models in the Z n circular lattice and square lattice in phase space. The operatorŜ i ≡ mŝ im is the total spin operator on the i-th site of phase space lattice since more than one spins can be confined in one lattice site. For the Z n circular lattice, the origin point (0, 0) in phase space is also a stable point, which can hold a central spin interacting with the spins on the Z n lattice. Therefore, the corresponding Heisenberg model up to a constant is whereŜ i(j) with i(j) = 1, 2, · · · , n represents the spins on the Z n lattice sites and S 0 is the central spin. This model is closely related to the central spin models [ 134,135]. The famous Mermin-Wagner theorem claims that a 1D or 2D isotropic Heisenberg model with finite-range exchange interaction can be neither ferromagnetic nor antiferromagnetic [136] at any nonzero temperature, which clearly excludes a variety of types of long-range ordering in low dimensions. Here in our model, the Flquet exchange interaction from point-like contact interaction has a Coulomb-like long-range behaviour. Therefore, the phase space lattices provide a possible platform to test the Mermin-Wagner theorem and search for other new phenomena with longrange interactions such as causality and quantum criticality [137][138][139][140][141][142][143], nonlocal order [126,144,145], etc.

Disorder & Localisation
By introducing a random temporal perturbation satisfying H (t + nT d ) = H (t) with T d = 2π/ω d to match the n : 1 nonlinear resonance [31,107], a random on-site potential in the tight-binding model can be produced, can be modulated by properly choosing the fluctuating function in H (t). For the Z n lattice, Eq. (58) is the 1D Anderson model on a ring. In the 1D Anderson model, the localization always happens no matter how weak the disorder is. In the phase space picture, the localised state is the superposition of the wave-packets around the Z n lattice site and with probability exponentially suppressed as the phase space distance increases. In the laboratory frame, the interpretation is that the probability density at a fixed position in the configuration space is localised exponentially around a certain moment of time [107], which is coined as the Anderson localisation in time domain.
A different but related scheme without utilizing classical nonlinear resonances is presented in [146] by Sacha and Delande. In this case the authors considered the problem of a single particle moving on a ring geometry under the action of a periodic force which however fluctuates randomly in time. The model is described by the following classical Hamiltonian, where the position of the particle on the ring is denoted by the angle variable θ and its canonical momentum by p. The Anderson problem is considered by introducing random disorder in space whereas here it is in time. Specifically, the time-dependent force consists of independent random Fourier components, e. g., with f k = f * −k independent random variables and ω d the driving frequency. Note that the spatial potential g(θ) is not necessarily disordered but any regular function. Inside the classical invariant torus one can employ the secular approximation and eliminate the fast oscillating part of the motion and arrive at an effective Hamiltonian describing the slow motion, Here P and Θ are canonical variables of the slow motion and they are related to the original ones by the transformation Θ = θ − ω d t and P = p − α − ω d . Assume the disordered potential satisfies the normal distribution |g k f −k | ∝ e −k 2 /(2k 2 0 ) with correlation length √ 2/k 0 and Arg(f −k ) are uniform random numbers in the interval [0, 2π). With a suitable choice of V 0 and k 0 , the eigenstates of (60) are Anderson localized even for energies higher than the standard deviation of the disordered potential. The Anderson localization in time crystals is revealed by measuring the probability that the particle passes a fixed position, which is localized exponentially around a certain moment of time. A possible experimental scheme for realising such Anderson localization is an electron in a Rydberg atom perturbed by a fluctuating microwave field [147]. The Anderson model was later extended to higher dimensions by considering particles moving in a three dimensional torus [148].
It is quite natural to combine interactions between particles with at the same time added temporal disorder, which extends the analysis of the Anderson localisation to interacting particle systems, i.e., the many-body localisation (MBL). In Ref. [131], Sacha et al proposed the following Bose-Hubbard Hamiltonian with on-site disorder where the interaction U ij and on-site disorder i are given by Eq. (59) and Eq. (55) respectively. The Hamiltonian (61) is valid provided the interaction energy N U ij (N is a total number of bosons) and the disorder j are much smaller than the energy gap. Although the disorder breaks the translational invariance of the phase space lattice, the system still possesses Floquet eigenstates with longer period 2πn/ω d equal to the period of disorder, i.e., H (t + nT d ) = H (t). The MBL in disordered Hubbard model X P Figure 11. Phase space lattice waves. The atoms (red dots) rotate around their equilibrium points unidirectionally, i.e., clockwise or anticlockwise, depending on the extremal properties (maximum or minimum) of the equilibrium points. The collective motions may arise from the phase space interaction between particles, which are the lattice waves in phase space.
has been observed in the experiment [149]. It is believed that MBL is constrained to systems with short-range interactions [150] but the long-range interaction alone does not exclude MBL [151]. The existence of MBL with long-range interactions is an important but unexplored problem. In the experiments with cold atoms, the interaction is short range since the atoms are neutral in general. The many-body Hamiltonian (61) from subresonantly driven protocol may present a new approach to solve this problem, as the long-range interaction U ij can be generated from contact interaction directly.
By modulating the strength of contact interaction with a proper temporal disorder, the particles can have an effective disordered potential in phase space. In such way, a two-particle molecule bound by Anderson localisation is created from the periodic driving [84].

Summary and Outlook
We started from Wilczek's time crystals breaking the CTTS of a static Hamiltonian [1,2], which has been disproved to exist in equilibrium systems [5,8,9,16,18]. However, in periodically driven (Floquet) systems, the time crystal behaviour does exist by breaking the DTTS [20][21][22][23][24], namely, the system responds at a fraction ω d /n of the original driving frequency. Such Floquet time crystal is a robust subharmonic mode oscillating in the laboratory frame with time period nT d , which is n times the driving period T d = 2π/ω d . Most works so far on time crystals are interested in stabilising the subharmonic mode via various mechanisms like MBL in spin chain models. This review paper, however, focuses on the many-body dynamics of subharmonic modes and the analogous condensed matter phenomena in Floquet time crystals [84,85,94,102,107,131,148]. Throughout this review, we choose the phase A gas of Floquet time crystals with spin chains. (Left) Four periodically driven spin chains (arrows with different colours) couple to an auxiliary spin (black arrow) via their rightmost spins. The circular arrows on the leftmost sides represent the periodic driving field and the noisy curves threading the spin chains represent the disorder to stabilise the subharmonic modes. (Right) Any two subharmonic modes of spin chains (coloured rounds) have the same effective interaction form U ij , which only depends on their oscillating phase φ i and φ j . The ensemble of subharmonic modes form a gas of Floquet time crystals.
space lattice picture to view the Flqouet time crystals or the subharmonic modes. Namely, we go to the rotating frame at the frequency of subharmonic modes, and place all the subharmonic modes in the phase space of the rotating frame according to their oscillating amplitudes and phases (or the action and angle variables). The resulting arrangement shows a periodic crystal structure in phase space, e.g., a circular Z n lattice and a 2D square lattice. The phase space lattice picture in the rotating frame [85,94,102], which is equivalent to the time crystal picture in the laboratory frame [84,107], is convenient for developing an analogous condensed matter theory as in solid state physics. We reviewed the single-particle band theory for phase space lattice [85,94] and the topological phenomena due to the noncommutative geometry of phase space [102]. We also discussed the effective interaction between two particles in phase space [102,127] beyond the single-particle picture. At this setting, the manybody Hamiltonian like the Hubbard model [84,107] and the Heisenberg model [102] can be realised in phase space lattice. By adding temporal disorder to the system, the Anderson model [107,[146][147][148] and the disordered Bose-Hubbard model [84,131] are also possible to be created in time crystals.
It is difficult to predict the research directions in the future although there are many unexplored topics in this new research field. Except mimicking the existing models in condensed matter physics, it would be more interesting to find something new in phase space lattice. For example, a phase space crystal should exist if the phase space lattice is strong enough to confine the atoms near its stable points. In this case, it is interesting to study the collective motion of atoms near their equilibrium points, which may form phase space lattice waves as indicated in Fig. 11. The lattice waves in phase space differ from the lattice waves in real space at a fundamental level that they are waves in noncommutative space. Each atom actually only rotates around its equilibrium position unidirectionally, i.e., either clockwise or anticlockwise depending on the extremal properties (minimum or maximum) of equilibrium points. Therefore, we conjecture the phase space lattice waves could have some topological properties, i.e., the nontrivial edge states [152].
Until now, the condensed matter phenomena in time crystals are only studied with the models of periodically driven quantum oscillators (i.e., driven quantum gas). One may ask if it is possible to develop a similar theory with periodically driven spin chain models. However, one should notice there is an important difference between these two models, i.e, each subharmonic mode of the driven oscillator can be a single-particle state while the subharmonic mode in the driven spin chain model is already a many-body state of spins. Nevertheless, we can investigate the manybody dynamics of subharmonic modes with spin models. As illustrated in Fig. 12, imagine there are several periodically driven spin chains interacting with an auxiliary centre spin on an equal footing. The subharmonic modes oscillate with frequency in the form of cos(ω d t/n + φ i ), where φ i is the oscillating phase of each subharmonic mode. The effective interaction between any two subharmonic modes should have the same form depending only on their oscillating phases. Then, the many-body dynamics of subharmonic modes can be described by the general Hamiltonian (46). In this way, the ensemble of subharmonic modes can be cast as a gas of Floquet time crystals. Some novel phenomena may arise in this model. For example, if the effective interaction is negligible, the time symmetry breaking processes of different spin chains are independent. However, if the interaction with centre spin is strong enough, the time symmetry breaking processes of different spin chains may affect each other and form a crystal in the oscillating phase φ i of time crystals, namely, from a gas-like state to a solid-like state.
The effects of the non-RWA Hamiltonian are also not fully understood yet. The phase space lattice is an emergent property of the RWA Hamiltonian, which is the leading Floquet-Magnus term in the rotating frame H RW A = H (0) F . The higher-order Floquet-Magnus terms from the non-RWA Hamiltonian [e.g., Eq. (38)] deteriorates the discrete lattice symmetry of the RWA Hamiltonian [e.g., Eq. (37)]. One can ask how the quasienergy band structure changes in the presence of non-RWA terms. In Ref. [85], it has been found that the band structure of the Z n lattice (31) is very robust for a large region beyond the RWA regime, but the tunnelling rate J l shows some sharp peaks at some detuning parameter, which has been identified as the resonant transitions induced by the fast non-RWA oscillating terms disregarded in the RWA [153]. In general, the current approach only applies for the stroboscopic dynamics of Floquet systems. The study on the non-RWA Hamiltonian describing the micromotions of system [46] is an important research direction in the future. Another interesting and fundamental research direction is the classification of possible spacetime groups for time crystals, which is analogous to the classification of space groups for static crystals. In Ref. [154], a complete classification of the 13 space-time groups in one-plus-one dimension (1+1D) is performed, and 275 spacetime groups are classified in 2+1D. Such study serves as the starting point for future research on the topological properties of time crystals.