Floquet time crystals in driven spin systems with all-to-all $p$-body interactions

A periodically driven Hamiltonian system exhibits discrete time-translation symmetry, that is, invariance of the Hamiltonian at time intervals separated by one period of the drive, $T$, $\hat{H}(t+T)=\hat{H}(t)$. A FTC is an out-of-equilibrium phase of matter emerging as a consequence of a physical observable breaking the discrete time translation symmetry of a many-body Hamiltonian [7, 6]. More precisely, the FTC phase can be defined by considering a class of initial states $\{|\psi\rangle\}$ and a generic choice of observable $\hat{O}$, such that the time-dependent expectation value in the limit of large system size $N$, given by

satisfies the following conditions: [9, 27]

1. 1.

   Time-translation symmetry breaking: $f_{O}(t+T)\neq f_{O}(t)$ while $\hat{H}(t+T)=\hat{H}(t)$

2. 2.

   Rigidity: $f_{O}(t)$ has a fixed oscillation period, without needing to fine-tune parameters in $\hat{H}$.

3. 3.

   Persistence: the oscillations of $f_{O}(t)$ persist for an infinitely long time

The conditions (1)-(3) imply that an FTC phase is not possible for generic chaotic systems, in which diffusion and equilibration typically preclude the existence of long-lived oscillations. Instead, FTC phases are expected to exist either in disordered systems (where diffusion can be suppressed by many-body localization [7, 8] or by other means [13, 28]) or in certain clean, integrable systems with highly regular motion. Clean Floquet time crystals have been studied in [27] for a Bose-Hubbard ladder, in [9] for the Lipkin-Meshkov-Glick model, and in [18] for a family of clock models. See also [29] for a precursor of these studies.

The emergence of an FTC implies an structural change in the nature of the Floquet states, which are the eigenstates of the unitary evolution operator $\hat{U}_{F}$ corresponding to $\hat{H}$ after one driving period. A key signature of this structural change is the emergence of eigenstate ordering [6], that is, the Floquet states become cat-like states, i.e., superpositions of macroscopically distinct states [7]. This is reflected in the reorganization of the Floquet spectrum, where an extensive portion of it is composed of groups of Floquet phases with spacing on the unit circle by an angle of $2\pi/q$, with $q$ and integer equal to the number of states entering in the superposition. As a consequence, the Floquet spectrum of $\hat{U}_{F}^{q}$ will be composed of clusters of $q$-fold degenerated Floquet phases, as all the Floquet phases in one of the groups seen in the spectrum of $\hat{U}_{F}$ will collapse to the same value. Eigenstate ordering has been recently proposed as the fundamental characteristic of an FTC phase [6], with the typical subharmonic response as its dynamical manifestation.

In the following we focus on a class of interacting magnetic systems called $p-$spin models [30, 31, 32]. Consider a system of $N$ spin-1/2 particles with $p$-body, $p\ll N$, Ising-like interactions, and in presence of an external homogeneous magnetic field, described by a time-independent Hamiltonian $H_{p}$. We will consider interactions to be homogeneous and of infinite range, i.e, all-to-all, as to avoid heating to a featureless “thermal” phase when the external driving is included ¹¹1In presence of short-range interactions, the full Hilbert space of the $N$ spin-$1/2$ particles is considered, and typically one requires the addition of disorder to induce localization and salvage some region of parameter space from the deleterious effects of heating [7, 8, 82]. The interaction Hamiltonian is permutationally symmetric and given an initial state in the symmetric subspace of the $N$ spin-$1/2$ particles, any time evolution will be restricted to that space. The Hamiltonian of the resulting family of $p$-spin models is

where $h$ is the strength of the external magnetic field, $\Lambda$ the strength of the $p$-body interaction, and we have introduced the collective spin operators $\hat{\mathbf{S}}=\sum_{i=1}^{N}\hat{\bm{\sigma}}^{(i)}/2$ with $\hat{\bm{\sigma}}$ a vector of Pauli operators. In terms of the collective spin operators, the symmetric subspace is spanned by the Dicke states, $\lvert S,M\rangle$, where $S=N/2$, resulting in a Hilbert space of dimension $N+1$, growing only linearly with the system size.

This family of models exhibits a rich and extensively studied variety of critical phenomena [30], including ground state quantum phase transition (GSQPT) between a paramagnetic and a ferromagnetic phase [34], excited state quantum phase transitions (ESQPT) ²²2For all the models with $p>2$, the phenomenology of the ESQPT is similar to that exhibited by the Lipkin model, see for example [24] as well as dynamical quantum phase transitions (DQPT) [36, 37, 38]. The properties of these phenomena are markedly different depending on the value of $p$, and in the following we summarize the aspects which are of particular importance to our construction of FTCs.

Critical phenomena in these models can be studied using a mean-field picture, which is exact in the thermodynamic limit $N\to\infty$ and describes the dynamics of the normalized mean spin $\mathbf{X}=\frac{\langle\hat{\mathbf{S}}\rangle}{S}$ on a unit sphere phase space. One tool to characterize this limit is given by the classical flow in phase space, $\dot{\mathbf{X}}=\mathcal{F}[\mathbf{X};h,\Lambda,p]$, which can be obtained via the Heisenberg equations motion for the spin operators after approximating $\langle AB\rangle\simeq\langle A\rangle\langle B\rangle$. Explicit expressions for this equations, which were previously studied in [37], are given in Appendix A. A second complementary tool is given by the semiclassical energy function $E(\mathbf{X};h,\Lambda,p)=\frac{\langle H_{p}\rangle}{S}$, which serves the role of a free energy, where the average is taken over a spin coherent state $|\theta,\varphi\rangle=e^{-i\varphi\hat{S}_{z}}e^{-i\theta\hat{S}_{y}}|S,S\rangle$. For the $p$-spin models this energy function takes the form

The critical points of different phase transitions correspond with certain structural changes of Eq. (3) as a function of the ratio $\frac{h}{\Lambda}$. In the limit $\frac{h}{\Lambda}\to\infty$ the semiclassical energy is a single well with a minimum at $Z=0$, indicating paramagnetic ordering. As the value of this ratio is reduced, Eq. (3) might undergo a saddle-node bifurcation at which a saddle and a local minimum (metastable ferromagnetic phase) emerge. This point is typically referred to as the spinodal point [39]. Upon further reduction of the value of $\frac{h}{\Lambda}\to 0$, the local minimum reaches a point of degeneracy with the current global minimum, signaling the GSQPT critical point [39], i.e., the point at which the ground state of the system changes character to ferromagnetic. The above scenario describes a first order GSQPT, as it is the case in the models with $p>2$ [30]. For the model with $p=2$, where one recovers the the Curie-Weiss paramagnet (a special case of the Lipkin-Meshkov-Glick or LMG model), the spinodal point and the GSQPT critical point coincide, and the transition from paramagnetic to ferromagnetic orderings is continuous or second order. For general $p$, the spinodal point can be shown to be (see Appendix B)

and the GSQPT critical point is given by

The behavior of the spinodal and ground state critical points is illustrated in the $p$-spin phase diagrams of Fig. 1a.

The emergence of a saddle point at the spinodal point implies the existence of a separatrix line in the classical phase space, which marks the boundary between two regions of distinct macroscopic motion of the mean spin. As a consequence of this, an extensive portion of the quantum spectrum is composed of states localized in the interior of the region enclosed by the separatrix line. In top-left panel of Fig. 2 we show an example of the phase portrait for the $p=2$-spin. The dark lines indicate trajectories on the inside of the separatrix and the light lines indicate trajectories on the outside of it. The set of eigenstates which are localized in the regions inside the separatrix correspond to the extensive portion of the spectrum which have undergone a clustering ESQPT [40, 41, 42, 24]. These localized states will be used to define FTC phases in the present work, since they will lead to breaking the discrete time translation symmetry of the driven $p$-spin system. Thus, we will consider exclusively the regime of $hW_{\rm GS}(p)<\Lambda$, where the system is in the ferromagnetic phase. Finally, we point out that $p$-spin models have been studied in the context of so-called Boundary Time Crystals [43, 44, 45], in which the the continuous time-translation symmetry is broken by the presence of dissipation. In our work, we focus exclusively on studying Floquet Time Crystals in these models.

Given an initial state localized inside the region enclosed by the separatrix of a $p$-spin model, we consider an external periodic drive in the form of a train of short pulses or “kicks”, as illustrated in Fig. 1b. The period of this train of pulses is given by $T$, which we will fix to $T=1$, and at every period the system undergoes an instantaneous rotation around the $x$-axis by an angle $\alpha_{\rm B}$. A single evolution step is thus given by

where we have taken $\hbar=1$. From the point of view of the quantum (finite-size) system, the external drive transforms the Hamiltonian system into a Floquet system. Conversely, in the mean-field description, the evolution of the mean spin becomes stroboscopic and dictated by an area-preserving map (APM)  [46] rather than a continuous phase space flow ³³3In fact, the mean spin evolves according to a “coordinate transformation” of the form $\mathbf{X}_{l+1}=\mathcal{A}[\mathbf{X}_{l};\alpha_{\rm B},h,\Lambda,p]$ with $l$ a discrete index accounting for the time step, the Jacobian being identical to $1$, hence the denomination of area preserving.. This latter fact is of utter importance, as the mean-field description in terms of an APM brings in new phenomenology [46] which is absent in the continous flow case. In particular, we will describe how the emergence of subharmonic system responses as a function of the drive parameters can be seen in phase space as the emergence of a resonance of the APM [48]. A resonance is a region of phase space enclosed by separatrices connecting hyperbolic periodic points. Trajectories in the interior of the resonance correspond to macroscopic motion of the mean spin in the form of oscillations exhibiting a strong periodic subharmonic component. See Fig. 2 for illustrations of this phenomenon. It is our aim to exploit the resonances of APMs in order to define FTC phases, and we will elaborate on this in the next section.

Finally, let us point out an important connection between the Floquet operator in Eq. (6) and the family of kicked $p$-spin models, recently introduced in Ref. [49]. In the case $h=0$, corresponding to a pure ferromagnetic $p$-spin, Eq. (6) exactly recovers the kicked $p$-spin unitary. More generally, the connection can be made naturally as long as $hW_{\rm GS}(p)\ll\Lambda$. This is precisely the regime in which the initial states discussed in Sec. II.2 are relevant. In this limit, we can rewrite the Floquet operator in Eq. (6) as

which has the form of a kicked $p$-spin model. The equality holds up to order $\mathcal{O}(h\Lambda)$ provided $\Lambda$ is not too large, and we defined the modified precession angle $\alpha(h)=\alpha_{\rm B}+h$. The latter will allow us to examine the robustness of the phase, for fixed value of $\Lambda$, using $h$ as tunable parameter ⁴⁴4We allow us to consider both positive and negative values of $h$, provided $|h|W_{\rm GS}\ll\Lambda$ holds. This change in sign changes the character of the $p$-spin ground state from ferromagnetic to anti-ferromagnetic, however the overall structure of the spectrum, i.e the different ESQPT lines remain invariant, as this change in sign can be seen as the application of $e^{i\pi\hat{S}_{y}}$. Thus, we still have an extensive portion of excited states undergoing a clustering ESQPT and defining a class of states which is localized at the interior of the region enclosed by the separatrix.. The APM for the mean-field limit of Eq. (7) has the form $\mathbf{X}_{l+1}=\mathcal{A}[\mathbf{X}_{l};\alpha(h),\Lambda,p]$, and explicit expressions are given in Appendix  A (see also [49] for details).

A resonance of an APM is a region of finite area in phase space delimited by separatrices connecting hyperbolic (periodic) points. Every resonance has a central periodic orbit with period equal to that of the resonance or some integer multiple of it [48]. Therefore, the trajectories in the interior of the resonance exhibit a strong subharmonic behavior with period $q$, with $q>1$ and integer. Consider a $1$:$q$ resonance, for an initial condition on or near the separatrix. Then, the stable or unstable character of the separatrix branches will only be apparent after $q$-steps of the stroboscopic evolution. Furthermore any initial condition inside the resonance will only look localized, as the parent trajectories of the classical flow of the Hamiltonian system, every $q$-steps of the stroboscopic evolution. For instance, take the system with $p=2$, the phase portrait for $\hat{H}_{p}$ is shown in upper-left panel of Fig. 2 (left sphere). In the undriven case, the motion of a state initially localized inside one of the lobes remains confined to that lobe. Once driving is included such that a resonance condition is met, the same type of initial state will now visit all other lobes before coming back to the initial one. In this case, $\alpha_{\rm B}=\pi$ and this process happens after $q=2$ steps. Similar scenarios occur for other values of $q$, see Fig. 2.

From this description of the classical map, we see that the macroscopic motion of our chosen set of initial states will be dominated by the central periodic orbit, and hence generic physical observables will display discrete time translation symmetry breaking. The existence of a resonance of the APM implies that the Floquet states, which in the limit $S\to\infty$ will correspond to trajectories inside the resonance, must have support on all lobes of the resonance, and hence could be written as superpositions of macroscopically distinct states, each of which has support inside one of the lobes of the resonance. Thus, a subset of Floquet states will look like ”cat-states”, and thus display eigenstate ordering in phase space, a key requirement needed for the definition of an FTC phase [6] ⁵⁵5If one finds Floquet states inside the resonance whose support is not in all lobes of the resonance, they will have a nontrivial evolution every step and thus will not satisfy an eigenvalue equation for the Floquet operator..

Using the resonances of an APM to define FTC phases brings in the robustness of the phase for free. That is, the observed subharmonic response will be robust to small changes in the system parameters, a statement that follows from the persistence of the fixed points. The latter concept can be explained as follows: Let us consider the APM $\mathbf{X}_{l+1}=\mathcal{A}[\mathbf{X}_{l};\alpha(h),\Lambda,p]$ giving the evolution of the mean spin in the thermodynamic limit. The periodic orbits on a resonance, hyperbolic or elliptic, are fixed points of $\mathcal{A}^{q}[\mathbf{X}_{l};\alpha(h),\Lambda,p]$ with $q$ the period of the resonance. A result in the theory of APMs guarantees the persistence of fixed points [46], that is, the resonance and its periodic orbits persist under the action of small perturbations, at most being moved to some new location in phase space, unless the tangent map evaluated at the fixed point has an eigenvalue equal to $1$. If the latter condition holds, the map undergoes a bifurcation and the periodic orbits might disappear.

It is desirable to have conditions on the system parameters for which emergence of a resonance is guaranteed. In the following we discuss two of these situations: resonance conditions and bifurcations in APMs.

Given an APM which describes the dynamics of the system, we can always assign a time-periodic Hamiltonian, $\hat{H}(t)=a\hat{H}_{1}+b\hat{H}_{2}\sum_{n=-\infty}^{\infty}\delta(t-nT)$, giving rise to the same APM equations in the thermodynamic limit. Then, a resonance condition is given by the values of $a,b$ for which the system receives an integer number of kicks in the form of $\hat{H}_{2}$ during the time it takes to complete a full cycle of the dynamics generated by $\hat{H}_{1}$.

For the Floquet system in Eq. (6) a resonance condition is given by the value of the angle of the drive for which the mean spin receives an integer number of kicks during the time it takes to complete a full period of the trajectory associated with the corresponding phase space flow. Correspondingly, for the kicked $p$-spin in Eq. (7), it is given by the situation at which the mean spin receives an integer number of kicks during the time it takes to complete a single precession, that is, when $\alpha_{\rm B}=\frac{2\pi}{q}$ with $q$ an integer. At a resonance condition, what used to be separatrices emerging from saddle points of the associated continuous flow, and enclosing regions of phase space where eigenstates are localized ⁶⁶6In the case of even $p$-spin models this localized eigenstates are symmetry broken states, breaking the $\mathbb{Z}_{2}$ symmetry of the Hamiltonian., now become separatrices connecting hyperbolic periodic points. The region enclosed by the separatrices becomes a $1$:$q$ resonance of the APM [48].

To explore the consequences of a resonance condition in our driven system we define the time-dependent Hamiltonian

where we have fixed the period of the drive to be $T=1$. Notice that this corresponds to the kicked $p$-spin Hamiltonian [49]. Eq. (8) can be brought into the form $\hat{H}(t)=\hat{H}_{\rm reso}+\hat{V}(t)$ (see Appendix C for details), where $\hat{H}_{\rm reso}$ is the resonance Hamiltonian and $\hat{V}(t)$ a time dependent perturbation which is typically of high frequency and whose effect vanishes on average.

The resonance Hamiltonian encodes the effects that the resonance condition induces in the Floquet system. Its form its derived in Appendix C, and it reads

where $\hat{\mathbf{O}}_{YZ}=(\hat{S}_{y},\hat{S}_{z})$ is the projection of the collective spin onto the $y$-$z$ plane, and $\vec{e}_{j}=\left(-\sin\left(\frac{2\pi}{q}j\right),\cos\left(\frac{2\pi}{q}j\right)\right)$ are the vertices of a $q$-regular polygon (see for instance [53] and chapter 6 of [54]). The resonance Hamiltonian is a sum of $q$ terms in the form of a $p$-twist each along one of the directions $\vec{e}_{j}$. The form of $\hat{H}_{\mathrm{reso}}$ implies that a region of finite size in the vicinity of the great circle in the $y$-$z$ plane develops a $1$:$q$ resonance, whose central periodic orbit has its points on the vertices of a $q$-regular polygon. This resonance then satisfies an emergent $q$-fold, $\mathbb{Z}_{q}$, rotational symmetry, which can be appreciated in the phase space portraits of Fig. 2 ⁷⁷7In fact the resonance has a symmetry dictated by the dihedral group $\mathrm{D}_{q}$. Extensive work investigating emergent symmetries at resonance conditions in some area preserving maps e.g delta-kicked harmonic oscillator, and their associated phase space structures, including crystalline and quasi-crystalline lattices, has been done by Zaslavsky, Sagdeev, Usikov and Chernikov [53, 54]..

One could, then, use every resonance condition of Eq. (7) to define an FTC phase, since an appropriate choice of the initial state and observable will reveal the subharmonic character of the system response, which will be robust to small changes in $h$ (leading to small changes in $\alpha_{\rm B}$ away from $\frac{2\pi}{q}$, as discussed in Sec. III.1). This is in fact the method employed in [23] and [19] for the identification of higher period FTC phases, which we have put into a more formal footing. However, the emergent symmetry discussed above is a feature of only a region of phase space in the vicinity of the $y-z$ plane. Following previous works [6], we require the emergent symmetry to be global in order to define a proper FTC. For mean-field models, this implies that the new symmetry affects the entirety of phase space (or a majority fraction of it). It is desirable, then, to introduce further requirements on the system parameters in order to position a given resonance condition as a point where a proper FTC phase exists. To this end we investigate the dynamics near the poles, i.e around the points $X=\pm 1$ and their bifurcations ⁸⁸8Notice that we do not need to analyze the bifurcations of a higher iteration of the map, as those higher period bifurcations will be seen, in the bare map, as points for which the eigenvalues of the tangent map are some root of unity, see for instance [46, 57]..

For an area preserving map, $1$:$q$ resonances can also emerge as consequences of a bifurcation process of a fixed point of the map. These, generally, are of the types described in the classification of generic bifurcations [57, 46]. Given a fixed point of the map, say $\mathbf{X}_{\rm fix}$, if the eigenvalues of $\left.\frac{\partial\mathcal{A}}{\partial\mathbf{X}_{l}}\right|_{\mathbf{X}_{\rm fix}}$ are equal to $1$, $-1$, or the $q$th root of unity, then the fixed point undergoes a tangent, period doubling, or period $q$ bifurcation, respectively. In the case of period doubling or period $q$ bifurcation, one also observes the emergence of a $1$:$q$, $q\geq 2$, resonance in phase space.

In [49] bifurcations of the fixed points at the poles $X=\pm 1$ were analyzed. It was shown that the kicked system with $p=2$ has a period doubling bifurcation at $\alpha_{\rm B}=\pi$, and all the other models with $p>2$ have $d$-$q$ bifurcations whenever $\alpha_{\rm B}=\frac{2\pi d}{q}$ with $d,q$ relative primes, $q\geq 2$ and $q>d$ (expressions for the tangent map eigenvalues are given in Appendix A). However, not all the allowed bifurcations in the classical limit introduced large structural changes with signatures in the quantum system. This can be understood by analysing the structure of the multi-body interaction term $\sim\hat{S}_{z}^{p}$. In terms of raising and lowering operators $\hat{S}_{\pm}=\hat{S}_{y}\pm i\hat{S}_{z}$, the interaction term will connect states in the $\hat{S}_{x}$ basis which are $l$-flips away, where $l\in\left[p,p-2,p-4,...,2(1)\right]$, and the set ends at $2(1)$ for even (odd) $p$-spins, that is, in the basis of $\hat{S}_{x}$ it introduces a coupling between states whose spin projection differ by $l$-units. Therefore, from all the discrete rotational symmetries emerging at resonance conditions of the term $\sim\alpha(h)\hat{S}_{x}$, only few of them are permitted by the interaction, which for the case of $p>2$ are also bifurcation points of the poles. From this considerations, we can define a reduced set of bifurcations

We refer to the resonance conditions which are also bifurcation points and are contained in $\mathcal{B}_{\rm bifu}^{(p)}$ as resonant-bifurcations. This set was previously identified in [58], where, using unitary perturbation theory, it was shown that large structural changes take place around these points ⁹⁹9Importantly, for the emergence of FTC phases one should also account for bifurcation points having multiplicities larger than one. If $\frac{2\pi}{m}\in\mathcal{B}^{(p)}_{\rm bifu}$, and $r\frac{2\pi}{m}<\pi$ for $r\in\mathbb{N}$, then that bifurcation point should also be included in $\mathcal{B}^{(p)}_{\rm bifu}$. In the case of even $p$-spins this is a consequence of the inherited $\mathbb{Z}_{2}$ symmetry, since it forces all odd period bifurcations to be double, thus they will essentially look as a even period bifurcation with twice the period, see [49] for details. For instance, when $p=6$ we have $\mathcal{B}^{(6)}_{\rm bifu}=\{\frac{2\pi}{6},\frac{2\pi}{4},\frac{2\pi}{2}\}$, however $2\times\frac{2\pi}{6}=\frac{2\pi}{3}<\pi$, we should include $\frac{2\pi}{3}$ in $\mathcal{B}^{(6)}_{\rm bifu}$..

At resonant-bifurcation points, not only a finite region in the vicinity of the $y$-$z$ plane exhibits an emergent $q$-fold rotational symmetry, but also these regions extend all the way up to the poles, making it a global symmetry of phase space. When the resonance condition is no longer satisfied, $\alpha_{\rm B}\notin\mathcal{B}_{\rm bifu}^{(p)}$, the linear term $\sim h\hat{J}_{x}$ introduces a tilt which displaces the resonance on the $y$-$z$ plane to a new location, and the time-independent part of the Hamiltonian now reads $h\hat{S}_{x}+\hat{H}_{\rm reso}$. Under these conditions, a wider region of phase space contains trajectories with a strong period-$q$ component. These come from both resonances: those coming from the resonance condition and those emerging as a consequence of the bifurcation, and from all the trajectories in between them. Even though not all phase space exhibits a $q$-fold rotational symmetry, a large portion of it does, and so it is an approximate symmetry. This approximate symmetry provides ground for the definition of an FTC phase.

Let us now collect all the different ideas presented in the previous sections into a single summary of a classical picture of Floquet time crystals in $p$-spin models.

Consider an initial state supported inside the ferromagnetic region of phase space, i.e the region enclosed by the classical separatrix of the phase space flow, of a $p$-spin model. When the system is driven, this region is transformed into a resonance of the APM, where now separatrices connect hyperbolic periodic points. The ensuing dynamics of the initial state is no longer confined inside the original support, but rather visits all the different lobes of the resonance. In fact, the system will undergo oscillations with the period equal to that of the central orbit of the resonance, returning to the region of the original support every $q$ applications of the drive, with $q$ the period of the resonance. Thus, the period of the resonance dictates the type of FTC behavior emerging as a consequence of the drive.

We have argued, and will provide evidence, that if the system has $p$-body interactions, with $p$ large enough, then by tuning the angle of the external drive to one of the values in the set $\mathcal{B}_{\rm bifu}^{(p)}$, FTC phases beyond a period double $2T$-FTC are accessible. These special angles are resonance conditions at which an emergent dynamical $q$-fold rotational symmetry arises in the system, leading to a special reconfiguration of Floquet states and Floquet phases, as we will discuss in next section. They are also bifurcation points, therefore the resonances which exist in phase space can only disappear via a second bifurcation process, providing a certain degree of robustness to the emerging FTC phases. The resonant-bifurcation points are then the key element we exploit in order to define $qT$-FTC phases, where $0<q\leq p$ as discussed above. In the remaining of the manuscript we will characterize, both for a finite system sizes and in the thermodynamic limit, the different $qT$-FTC phases emerging in periodically driven $p$-spin models.

Eigenstate ordering of the Floquet states refers to a special structure of the eigenvectors of the Floquet operator $\hat{U}_{F}$, which have the form of cat-like states, i.e superpositions of macroscopically distinct states [7, 6]. In the mean-field limit, this translates to Floquet states which have support in all of the $q$ lobes of a $1$:$q$ resonance, that is, superpositions of $q$ states each of which is localized in one of the lobes of the resonance. Therefore a portion of the Floquet spectrum is composed of groups of $q$ Floquet phases, where each of the Floquet phases in one of these groups define the vertices of a $q$-regular polygon in the unit circle.

In the case of $p=2$ and $\alpha(h)\simeq\pi$, i.e the FTC phase in the LMG model [9], this phenomenon was referred to as $\pi$-pairing of Floquet phases. More precisely, the Floquet states with support inside the $1$:$2$ resonance come in pairs, with each member of a pair having a definite parity under $e^{i\pi\hat{S}_{x}}$. The eigenphases of Floquet states in one of these pairs are located at diametrically opposed points on the unit circle and thus they differ by $\pi$. As a consequence, when considering the operator $\hat{U}_{F}^{2}$ the Floquet phases corresponding to Floquet states in those pairs will be degenerate. Thus, the spectral statistics of $\hat{U}_{F}^{2}$ will deviate from the expected Poissonian behavior ¹⁰¹⁰10Here Poisson statistics of the adjacent level ratios is expected as we are in a parameter regime where classical system is regular. due to a strong degeneracy arising from to the clustering of levels by pairs.

More generally in FTC phases of period $q$ with $q>2$ we expect to see $\frac{2\pi}{q}$-pairing of the Floquet phases of $\hat{U}_{F}$ as the manifestation of eigenstate ordering. This pairing will lead to a strong degeneracy in the Floquet phases of $\hat{U}_{F}^{q}$, with its spectral statistics deviating from that of a Poisson distribution. In Fig. 3 we show an example of the Floquet phases for the system with $p=3$ around the region of the $3T$-FTC phase. Notice the emergent apparent regularity of the spectrum in Fig. 3a, and how the $\frac{2\pi}{3}$-pairing is revealed through the degeneracy of the Floquet phases of $\hat{U}^{3}_{F}$ in Fig. 3b. The existence of this degeneracy is a signature of an emergent dynamical symmetry, a discrete $q$-fold rotational symmetry, which is also manifested in the form of the effective Hamiltonian of $\hat{U}_{F}^{q}$, give by $h\hat{S}_{x}+\hat{H}_{\rm reso}$ as discussed in Sec. III.2. To describe this behavior we will use the average adjacent spacing ratio [61], a standard measure of spectral statistics, defined by

where $d_{j}=\mu_{j+1}-\mu_{j}$ is the eigenphase spacing and $N=2S$. Further, we define the normalized average adjacent spacing ratio $\tilde{r}=\overline{r}/r_{\rm POS}$ where $r_{\rm POS}$ is the value for a Poisson distribution, hence any value of $\tilde{r}<1$ will indicate a high degree of eigenphase clustering, signaling the onset of the FTC phase.

As stated above, the emergence of $qT$-FTC can be diagnosed by the clustering of an extensive portion of the eigenspectrum of $U_{F}^{q}$. This phenomenon admits an interesting connection to a clustering excited-state quantum phase transition (ESQPT) [24]. For a time-independent Hamiltonian with a one degree of freedom, a clustering ESQPT can be identified by a series of vanishing energy gaps in the spectrum, forming a continuum in the thermodynamic limit. Due to the closing gaps, the energy levels cluster around these regions, and the density of states diverges [40, 41, 24]. ESQPTs have been mostly identified in systems with clear-cut classical limits (like the $p$-spin models), and their existence is intimately connected to the emergence of separatrices in the mean-field phase space. Even though the $p$-spin models present ESQPTs, the phenomena shown here in the driven $p$-spin models is different, since the associated classical phase spaces are different. In general, the emergence of ESQPTs in driven quantum systems is a relatively unexplored phenomenon, with only few examples [62, 63, 64].

From a pure kinematic point of view, the nonFTC-FTC transition, as seen by the emergence of eigenstate ordering, displays all the features of a clustering ESQPT [41]. For a system whose thermodynamic limit is described by a single degree of freedom, a clustering ESQPT corresponds to a singularity of the density of states in the form of a logarithmic divergence. This divergence originates in the diverging oscillation periods of trajectories corresponding to initial conditions which get closer and closer to a separatrix line as a function of a control parameter. Therefore, when one examines the energy spectrum, it is crossed by a curve along which excited states cluster together towards degeneracy. For our system this observation follows naturally from the mean-field phase space analysis made in Sec. III, as the nonFTC-FTC transition, in this limit, is accommodated by the emergence of resonances of the APM. Let us assume a $1$:$q$ resonance leading to a $qT$-FTC, recall that under the effective representation of the dynamics every $q$ steps the partial separatrices connecting hyperbolic periodic points behave as actual separatrices emerging from saddle points. In fact, they correspond to separatrices of the effective Hamiltonian of $\hat{U}_{F}^{q}$, given by $\hat{H}_{\rm eff}(h)=h\hat{S}_{x}+\hat{H}_{\rm reso}$, with $h$ playing the role of control parameter. As such the density of states will display a logarithmic divergence, a phenomenon we discuss in more detail in Appendix  D.

The structural change of an extensive portion of the Floquet states from lacking eigenstate ordering to displaying eigenstate ordering has a direct consequence in the dynamics. In fact, the emergence of an extensive $\frac{2\pi}{q}$-pairing of the Floquet phases makes $\omega=\frac{2\pi}{q}$ to be the dominant frequency in the system response. As such, the macroscopic motion of the mean spin undergoes a drastic change, going from approximating Larmor precession at some frequency, whose value is in principle arbitrary and changes continuously with the system parameters, to be precessing with a frequency locked at the $q$th subharmonic of the frequency of the drive. This change in the macroscopic motion can be thought of as a dynamical quantum phase transition (DQPT) [65, 66, 26], and technically its associated with an order parameter DQPT, or DQPT of type-I [25]. We notice that, the drastic change in the macroscopic motion concerns only the frequency of precession and not the precession axis, therefore the asymptotic value of the quasi-steady state of an observable which begins out of equilibrium remains unchanged, having the consequence that the “usual” DQPT order parameters, such as the time-averaged magnetization and the time-averaged two point correlation function, fail at detecting this dynamical transition.

In order to have a faithful identification of this DQPT we resort to higher-order correlation functions, in particular out of time order correlators (OTOCs). OTOCs are widely used in studies of out-of-equilibrium dynamics and information scrambling [67, 68, 69]. Recently, it has been shown that in cases where a dynamical order parameter might not be well defined or it is difficult to construct, the long-time average of the OTOC can be used to dynamically detect quantum phase transitions, be it of equilibrium or dynamical nature [70, 71, 72]. The OTOC is defined as

where $\hat{W}(0)$ and $\hat{V}$ are two operators which commute and $\hat{W}(t=lT)=\hat{U}_{F}^{\dagger l}\hat{W}(0)\hat{U}_{F}^{l}$ is the Heisenberg evolution of $\hat{W}(0)$, and the average is taken on some reference initial state. We consider the infinite time average of Eq. (12):

We will be focusing on OTOCs with operators $\hat{W}(0)=\hat{V}=\frac{\hat{S}_{z}}{S}$, which we denote $\mathrm{F}_{Z,Z}(t)$ and $\mathrm{F}_{Z,Z}^{\infty}$, where the average is taken over an “infinite temperature” state $\rho_{0}=\frac{\mathbb{I}}{N+1}$. Using the long-time average of the OTOC as a diagnostic, we will see that the transition between non-FTC to FTC behavior can be regarded as a DQPT, where $\mathrm{F}_{Z,Z}^{\infty}=0$ in the non-FTC phase and $\mathrm{F}_{Z,Z}^{\infty}\neq 0$ inside the FTC phase.

To complement the characterization of the FTC phases made via the spectral statistics and the OTOC, we study phase diagrams of our FTC phases in the thermodynamic limit. These will be constructed based on the dominant frequency of the mean response of the system, obtained via the phase space average of the long time evolution of the classical $ZZ$ correlation function. This correlation function is defined as

where $Z_{l}$ is the stroboscopic evolution of the $z$-coordinate of the mean spin, and the brackets $\langle.\rangle_{\rm PS}$ indicate phase space average. Furthermore, to directly observe the subharmonic character of the system response we investigate the power spectrum $|\mathcal{C}^{Z}_{\omega}|^{2}$ of the averaged correlation function, where $\mathcal{C}^{Z}_{\omega}$ is given by

the discrete Fourier transform of the finite-time averaged correlation performed over $T_{\rm max}\gg 1$ periods. A phase diagram as a function of $\alpha(h)$ and $\Lambda$ can then be constructed by identifying the regions in parameter space where $|\mathcal{C}^{Z}_{\omega}|^{2}$ presents a peak at $\omega=2\pi/q$. The actual quantity to be computed reads:

where the normalization is included for convenience, since for systems with coexisting FTC phases, the $2T$-FTC might overshadow all the other phases.

As mentioned in the previous subsection, the change in the macroscopic motion only concerns the frequency of oscillation, as such the asymptotic value of classical two point correlation functions in the dynamical steady state, remains unchanged. However, we are constructing phase diagram based on the frequency of this correlation. Furthermore, the phase space average will, in absence of a single strong global frequency, acts as a classical dephasing process. However, inside an FTC we expect a large portion of phase space to be populated by trajectories whose dominant frequency is the $q$th subharmonic of the drive frequency. Therefore, the existence of an non-negligible peak at $\omega=\frac{2\pi}{q}$ in the phase-space-averaged system response provides strong evidence of the extensivity of the FTC under study.

The $2T$-FTC phase emerging at the $1$-$2$ bifurcation in the vicinity of $\alpha_{\rm B}=\pi$ is the most prominent one for all even values of $p$, as the emergent symmetry agrees with the $\mathbb{Z}_{2}$ symmetry the Floquet system inherits from the time independent Hamiltonian. In particular, the case with $p=2$ recovers the FTC phase in the LMG first discussed in [9].

In Fig. 4a we show the results of the characterization of the $2T$-FTC phase for two example cases with $p=2,4$. Fig. 4a.1.1 shows a snapshot of $f_{Z}(t)$, computed with the initial state $|\theta,\varphi\rangle=|\pi/5,0\rangle$ in the long time limit. The period doubled oscillation is clearly seen for both the cases of $p=2$ and $p=4$, indicated by circles and triangles respectively. Fig. 4a.1.2 shows the normalized power spectrum of $f_{Z}(t)$, displaying the clear peak at $\omega=\frac{2\pi}{2}$, accordingly to the subharmonic response of the system in the $2T$-FTC phase. The extension and robustness of these 2T-FTC phases are studied in the subsequent figures. Fig 4a.2 shows the normalized mean adjacent spacing ratio $\tilde{r}(\alpha(h),\Lambda)$ in the vicinity of $\alpha_{\rm B}=\pi$ for systems with $p=2,4$, respectively. The nonwhite areas indicate the parameter region where the spectral statistics deviates from that of a Poisson distribution due to large degeneracies by pairs in the spectrum of $\hat{U}_{F}^{2}$. Finally, Fig. 4a.3 shows the long-time average of the OTOC $\mathrm{F}^{\infty}_{Z,Z}(\Lambda,\alpha(h))$ for systems with $p=2,4$. The green area indicates $\mathrm{F}^{\infty}_{Z,Z}=0$ and thus non-FTC phase, whereas nongreen area indicates $\mathrm{F}^{\infty}_{Z,Z}\neq 0$ thus FTC phase.

At small values of $\Lambda$ the phase boundary of the FTC is essentially a straight line with a slope proportional to $\Lambda$, as can be seen in Fig. 4a.2. At larger values of $\Lambda$, the shape of this boundary becomes more complicated, however it can be computed using the APM in the mean-field limit. We give an example for the system with $p=2$ in Appendix  A. The extent of the $2T$-FTC phase in the direction of $\Lambda$ is, in principle, unconstrained and there will always be a narrow strip of FTC phase around $\alpha_{\rm B}$ in the limit $\Lambda\to\infty$. This can be understood from the fact that the Floquet system at $\alpha(h)=\alpha_{\rm B}$ is fully integrable and never transitions to chaos. In fact, the largest Lyapunov exponent in the limit of strong chaos (thermal regime), is given by $\lambda_{+}(\alpha(h),\Lambda,p)=\ln\left[\Lambda(p-1)\sin(\alpha(h))\right]-(p-1)$ (see [49] for details). Hence, as $\alpha(h)\to\pi$ one requires $\Lambda\to\infty$ in order to give $\lambda_{+}>0$. This fact is illustrated in Fig. 4a.2,a.3 (also Fig. 6) where the region of the FTC phase appears to extent indefinitely in the $\Lambda$ direction. This is a feature only of the $2T$-FTC phase and its emergence in the vicinity of $\alpha_{\rm B}=\pi$, and it is not present in the other FTCs present in this system.

As discussed in Sec. III, we can go beyond period doubling FTCs by setting $\alpha_{B}=\frac{2\pi}{3}$ leading to a $1$:$3$ resonance and identifying for which models a significant period-tripling bifurcations takes place, i.e. for which model $\frac{2\pi}{3}\in\mathcal{B}_{\rm bifus}^{(p)}$. The first of those models corresponds to $p=3$, the second and third ones have $p=5,6$, respectively. We illustrate the characteristic period tripling behavior of this phase using the systems with $p=3,6$ in Fig. 4b. In Fig. 4b.1.1 we show a snapshot of $f_{Z}(t)$ computed with the initial state $|\phi_{0}\rangle=|0,0\rangle$, where dots and triangles show $p=3,6$ respectively. The period tripling behavior is manifested in the normalized power spectrum of $f_{Z}(t)$ in Fig. 4b.1.2, where the peak is locked at $\omega=\frac{2\pi}{3}$ for both cases.

The extension of the phase in parameter space is investigated using the normalized average level spacing ratio $\tilde{r}(\Lambda,\alpha(h))$ and the long time average of the OTOC $\mathrm{F}^{\infty}_{Z,Z}(\Lambda,\alpha(h))$. In Fig. 4b.2 we show $\tilde{r}(\Lambda,\alpha(h))$ in the vicinity of the bifurcation point, where white areas indicates Poisson statistics and nonwhite areas indicates deviations from Poisson statistics, which we connect to presence of FTC behavior. Notice that the FTC phase for the system with $p=3$ in Fig. 4b.2.1 is highly constrained, as chaos emerges around this value of $\alpha(h)$ fairly rapidly  [49], erasing any footprint of the locked subharmonic periodic behavior. On the other hand, for the system with $p=6$ the phase is wider in both, the $\alpha$ and $\Lambda$ directions. For this case, the system gains stability thanks to the $\mathbb{Z}_{2}$ symmetry of the even $p$ spin systems, and the period-3 bifurcation is double, thus being structurally similar to a $1$-$6$ bifurcation. The FTC phase extents to larger values of $\Lambda$ since higher values of $p$ have a delayed emergence of global chaos [49]. The long time average of the OTOC in Fig. 4b.3 shows very good agreement with the results of $\tilde{r}$ for the $3T$-FTC phase. Here, green indicates nonFTC behavior and nongreen indicates FTC behavior.

Following the discussion in Sec. III, models where $\frac{2\pi}{4}\in\mathcal{B}_{\rm bifus}$ host a $4T$-FTC phase in the vicinity of $\alpha_{\rm B}=\frac{2\pi}{4}$. We fix the angle of the drive at this value and investigate the emergence of this phase with two exemplary systems with $p=4,6$. Fig 5a.1.1 shows, for the initial state $|\theta,\varphi\rangle=|\pi/5,0\rangle$, a snapshot of $f_{z}(t)$ in the long time limit, revealing the period quadrupling behavior (circles and triangles indicate $p=4,6$, respectively). Correspondingly, Fig. 5a.1.2 shows the normalized power spectrum of $f_{Z}(t)$, displaying a clear peak at the subharmonic frequency $\omega=\frac{2\pi}{4}$ in both cases.

The extension of this FTC phase in parameter space is investigated with both the normalized average level spacing ratio $\tilde{r}(\Lambda,\alpha(h))$ and the long-time average of the OTOC $\mathrm{F}^{\infty}_{Z,Z}(\Lambda,\alpha(h))$. The normalized spectral statistics indicator shows a region in parameter space associated with the $4T$-FTC phase, which grows in the $\Lambda$ direction as we increase the value of $p$, compare Fig. 5a.2.1 ($p=4$) and Fig. 5a.2.2 ($p=6$), here white indicates Poisson statistics and nonwhite indicate a value deviating from Poisson statistics which we associate with some degree of level clustering by quartets. The long time average of the OTOC is shown in Fig. 5a.3 for the systems with $p=4,6$, respectively. Here green indicates $\mathrm{F}^{\infty}_{Z,Z}=0$ and nongreen indicates $\mathrm{F}^{\infty}_{Z,Z}\neq 0$. Notice that $\mathrm{F}^{\infty}_{Z,Z}$ agrees with $\tilde{r}$ on the region they associate with signatures of a $4T$-FTC. We point out that these two indicators identify very well the finite size phase diagram of our FTC phases, regardless of the period of the phase.

An important issue arises here. Could we have identified a $4T-FTC$ for a system with $p=2$? As we mentioned in the Sec. I and Sec. III, area preserving maps naturally develop resonances, a simple way of identifying them is via the resonance conditions associated with a delta-kicked Hamiltonian. Thus, one could define FTC phases (or, rather, their precursors) using those resonances. However, the requirement of a global $\mathbb{Z}_{q}$ dynamical symmetry is usually not satisfied. This is the case for $p=2$ when one tries to define a $4T$-FTC at at the resonance condition $\alpha_{\rm B}=\frac{2\pi}{4}$. We discuss this case in more detail in Appendix E.

So far we have discussed FTC phases up to period quadrupling. However, higher period FTC phases are accessible provided we work with a system with large enough $p$-body interactions. A given kicked $p$-spin will always have a $1$-$p$ bifurcation at $\alpha_{\rm B}=\frac{2\pi}{p}$, whose vicinity can be used to define a $pT$-FTC phase. Additional $qT$-FTC phases with $q<p$ could be defined, in the vicinity of the appropriate value of $\alpha_{\rm B}$, if $\alpha_{\rm B}=\frac{2\pi}{q}\in\mathcal{B}_{\rm bifu}$.

Let us consider the system with $p=6$, whose lower period FTC phases have been analyzed in the previous subsections. We focus on values of $\alpha(h)$ in the vicinity of $\alpha_{\rm B}=2\pi/6$, point at which the corresponding mean-field system has a $1$-$6$ bifurcation. The response of the system is characterized using $f_{Z}(t)$ for an initial state of the form $|\theta,\varphi\rangle=|\pi/10,0\rangle$. Fig. 5b.1.1 shows a snapshot in the long time limit of this response, where a clear $6$-periodic oscillation can be observed. To complement this observation, we look at the normalized power spectrum of the $f_{Z}(t)$ signal in Fig. 5b.1.2, showing a clear peak at the frequency $\omega=2\pi/6$, signaling a $6T$-FTC.

To investigate the extension of this FTC phase in parameter space we look at both, the average adjacent level spacing ratio $\tilde{r}(\Lambda,\alpha(h))$ and the long-time average of the OTOC $\mathrm{F}^{\infty}_{Z,Z}(\Lambda,\alpha(h)$. In Fig. 5b.2.1 the spectral statistics indicator is shown, where nonwhite areas indicate deviation from Poisson statistics providing evidence for the level clustering by sextets, a signature of the onset of the $6T$-FTC phase. In Fig. 5b.3.1 the long time average of the OTOC is shown. Here, green indicates $\mathrm{F}^{\infty}_{Z,Z}=0$ and nongreen indicates $\mathrm{F}^{\infty}_{Z,Z}\neq 0$ signaling a nonFTC and FTC phases, respectively.

We close the analysis of finite-size systems by stressing that the two novel indicators proposed here, one based on the spectral statistics of the appropriate power of the Floquet operator and the other based on the behavior of the OTOC, have shown to provide excellent characterizations of the nonFTC-FTC phases. These two indicators are also complementary, as the first one is purely kinematic and the second one purely dynamic, providing then interesting tools for the characterization of FTC phases in systems beyond the $p$-spin models studied here.

The results in previous subsections characterize the emergent FTC phases for finite system sizes. In this subsection we will construct the phase diagram of the different FTC phases directly in the thermodynamic limit, for kicked $p$-spin models up to $p=6$. For these models, this can be done simply by solving numerically the nonlinear classical equations of motion for the magnetizations $(X,Y,Z)$. From this, we can compute the measure $\mathcal{G}(\Lambda,\alpha(h))$ introduced in Sec. IV.3. In the following, we show results where the phase space average was done using an uniform grid on the sphere with $14000$ points, and evolution times going up to $T_{\rm max}=16000$.

Results are shown in Fig. 6. In all cases with even $p$ the $2T$-FTC phase emerging around $\alpha_{\rm B}=\pi$ is robust in the $\Lambda$ direction. As pointed out previously, for this angle of the drive the system never transitions to chaos, and thus the system avoids thermalization in the long-time limit (see discussion in Sec. V.1). The $4T$-FTC phase in the system with $p=4$ exists up to $\Lambda\sim 5$ as it is then that the system becomes globally chaotic [49]. On the other hand, as we increase $p$ the $4T$-FTC phase gains in extension along the $\Lambda$ direction, since higher values of $p$ delay the emergence of global chaos in this family of models [49].

Some other important observations can be drawn from these results. For the models with odd values of $p$, FTC phases tend to be highly constrained (notice the extent of the horizontal axis in Figs. 6 (b) and (d) with respect to the rest). For instance, for $p=3$ the model hosts a $3T$-FTC phase around $\alpha_{\rm B}=\frac{2\pi}{3}$, however the poles at the bifurcation point are unstable and chaos emerges fairly rapidly for this model around this value of the angle of the drive. In fact as seen in Fig. 6b the phase only goes up to $\Lambda~{}\sim 1.5$. The $3T$-FTC phase in the model with $p=5$ is essentially identical to that in the model with $p=3$, whereas the $5T$-FTC is narrower in the direction of $\alpha$. This phenomenon keeps occurring as $p$ becomes larger in the models with odd $p$’s. Finally, the $3T$-FTC phase in the system with $p=6$ is structurally similar to the $6T$-FTC phase, see Fig. 6e. This is due to the fact that the APM of the kicked $p$-spin, $\mathcal{A}[\mathbf{X}_{l};\alpha(h),\Lambda,p]$, is a double reversible map for models with even $p$ [46, 49], in other words, the inherited $\mathbb{Z}_{2}$ symmetry imposes special constraints to bifurcations which have odd period, in this case they are double, meaning that one observes the emergence of two $1$:$3$ resonances looking structurally identical to a single $1$:$6$ resonance. Finally, we highlight the resemblance between the finite size phase diagrams constructed with the spectral statistics indicator and the long-time average of the OTOC, and shown in Fig. 4 and Fig. 5, with those constructed for the systems in the, $N\to\infty$, thermodynamic limit, discussed in the present subsection and shown in Fig. 6.