Time Crystals in the Driven Transverse Field Ising Model under Quasiperiodic Modulation

We investigate the transverse field Ising model subject to a two-step periodic driving protocol and quasiperiodic modulation of the Ising couplings. Analytical results on the phase boundaries associated with Majorana edge modes and numerical results on the localization of single-particle excitations are presented. The implication of a region with fully localized domain-wall-like excitations in the parameter space is eigenstate order and exact spectral pairing of Floquet eigenstates, based on which we conclude the existence of time crystals. We also examine various correlation functions of the time crystal phase numerically, in support of its existence.


I. INTRODUCTION
Our understanding of the out of equilibrium phase structures of periodically driven (Floquet) quantum many-body systems has made impressive progress over the past decade, see for example the review [1] and the references therein. Among the many interesting discoveries made in this context, the prediction and the consequent observation of a time-crystalline phase predicted by Wilczek [2][3][4][5][6][7] in Floquet systems has stimulated a lot of interest. The first concrete proposal to circumvent the no-go theorem by Oshikawa and Watanabe [8], assessing the impossibility to have spontaneous breaking in equilibrium systems, was to consider periodically driven systems. Floquet time-crystal were predicted theoretically in [9,10] and soon after experimentally observed in [11,12]. In the last few years the literature on timecrystals has grown enormously in many different directions, see for example [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29]. A recent overview can be found in Ref. [30].
Unlike previous works on manipulating and engineering effective Hamiltonians in the prethermal stage of Floquet evolution building upon high-frequency expansion [31] and on the edge states in Floquet systems relating to non-trivial topology of bulk bands, the long-time state of the discrete time crystals (DTC) is characterized by persistent oscillation of local order parameter, which implies the existence of non-trivial bulk spatio-temperal order. On the eigenstate level, it manifests eigenstate order in that almost all Floquet eigenstates possess nontrivial long range order.
Generic non-integrable Floquet systems will be heated up to the infinite temperature ensemble due to persistent pumping of energy by the driving [32][33][34], while the long time behavior of integrable Floquet systems is described by the periodic Gibbs ensemble [35][36][37]. This fact implies that the long time steady states of both cases always have trivial correlations, thus excludes any non-trivial phase * On leave structures. A way out, as first proposed in Ref. [9], is to consider many-body localized (MBL) systems [38,39]. MBL in the presence of high-frequency periodic drive was considered as well [40,41]. Many-body eigenstates may be classified by their broken symmetries and topology in the MBL phase [39,[42][43][44], analogous to that of phases transitions in equilibrium systems. Following this line, Ref. [45] showed that the driven quantum Ising chain hosts a Z 2 time crystal phase.
Localization can be induced not only by quenched disorder, but also by quasiperiodic modulation of couplings. As first discussed by Azbel [46], Aubry and Andre [47] and their generalizations [48][49][50][51][52], incommensurability leads to a localization-delocalization transition even for the 1D tight-binding lattice model. In this work, we study the out of equilibrium phase structure of a driven transverse field Ising model (TFIM) with quasiperiodic (QP) modulation of Ising coupling. We work with an integrable model which allows us to use various analytic tools and present reliable numerical results. Making use of the Jordan-Wigner (JW) transformation, the driven TFIM is mapped to a driven Majorana chain. We start by analyzing the Majorana edge modes and localization of single-particle excitations and show that localized domain-wall-like excitations lead to eigenstate order and time crystal order, therefore building a connection between localization of single-particle excitations and time crystals.
The paper is structured as follows. In Sec. II we introduce our model and discuss the implication of its symmetries on the phase structure. In Sec. III we present our results regarding the Majorana edge modes, localization of single-particle excitations, long-range order of excited states in order. Based on these we conclude the existence of time crystal phase and present numerical results supporting our statement. We summarize our results in Sec. IV.

II. MODEL
The driven QP-TFIM we consider consists of a twostep drive protocol, which is described by the following time-periodic Hamiltonian with period Here σ x,y,z j are the Pauli operators acting on site j and the Ising couplings are given by a smooth function, which we take of sinusoidal form Quasiperiodicity implies that the wavelength 1/Q is incommensurate to the lattice constant a, namely Q is an irrational number (we set a = 1 hereinafter). To analyze such a system theoretically, one can start by approximating the modulation by a sequence of periodic functions, which retrieve translational invariance. In this work we choose Q as the golden mean Q = ( √ 5+1)/2 and approximate it by the consecutive ratios of Fibonacci numbers q n+1 /q n , where the Fibonacci sequence is given by the recurrence relation q n+1 = q n + q n−1 and the initial conditions q 1 = 1, q 2 = 1.
The undriven version of this QP-TFIM, with the same form of Ising couplings J j , was studied in great detail in Ref. [53], where it was found that in the strong modulation regime J < A J new gapless phases with localized or multifractal excitations emerge. In a similar way, the constant Ising coupling J of our driven model controls the density of weak couplings and, as we will see below, plays a crucial role in determining the localization properties of (1). We have confirmed in our numerics that the phase ϕ J is irrelevant to magnetic order and localization properties of single-particle excitations away from the phase boundaries, similar to what was observed in the undriven case [53,54].

A. Symmetries of the Floquet operator
Our interest will be in the long-time behavior of the driven system, so it suffices to consider the Floquet operator of (1), which reads There are several symmetries of this Floquet operator. The Z 2 Ising symmetry, denoted by the symmetry operator P = j σ z j , is inherited from the Hamiltonian (1) and satisfies P U F P † = U F . Furthermore, the Floquet operator U F simply changes sign under the shift T J J → T J J + π or T b b → T b b + π, which can be compensated by shifting the quasienergy → + π while leaving the Floquet eigenstate unchanged. This implies that all the physical properties have a periodicity π in both T J J and T b b. Finally, we consider the reflections T b b → −T b b and T J J → −T J J. For the T b b reflection, U F preserves its form after applying local rotations to all the sites, such that σ z j → −σ z j . Instead, the T J J reflection can be combined with ϕ J → ϕ J + π and local rotations on the even (or odd) sites, such that σ x j → −σ x j . Taking into account all these symmetries, we can restrict our discussion to the rectangular region [0, π]×[0, π] in the T b b-T J J plane. Another simple observation is that we can assume A J > 0, as a change of sign is equivalent to a phase shift ϕ J → ϕ J + π.

B. Majorana representation
To describe the time evolution in terms of noninteracting fermions, we introduce a pair of Majorana fermions for each spin-1/2 The Majorana operators satisfy the anticommutation relation {γ i , γ j } = 2δ ij . Then, the Floquet operator can be expressed as Note also that there is a boundary term if periodic boundary conditions are imposed on the spin chain. However, this makes no trouble in our case as we will consider a semi-infinite chain in the study of edge modes, while for sufficient long chains the boundary term is marginal in detecting bulk properties. The Majorana representation is superior to the fermionic one in that it allows us to treat the problem of edge modes and localization of excitations on equal footing. A crucial feature of the Floquet operator is that it consists of two unitaries acting on disconnected Majorana pairs. We will see that this brings further simplifications, allowing us to obtain analytic expressions for some of the phase boundaries. For the moment, we derive the eigenequations for the Fermionic eigenmodes, which are defined in terms of the Majorana operators as and satisfy U † F Γ s U F = e −i s Γ s . Since shift of the quasienergy s → s + 2π results in the same eigenmode, we can restrict the values of s to the first Brillouin zone [−π, π]. In general the Γ s (Γ † s ) are fermionic annihilation (creation) operators, satisfying {Γ s , Γ † s } = δ s,s . Since Γ s , Γ † s are associated to opposite quasienergies, in the following we further restrict s ∈ [0, π]. The action of the two unitaries U 1,2 on a single Majorana fermion can be easily written down [see Eq. (A1)] and, with j in the bulk, yields Finally, based on Eq. (7) we discuss the extension of the weak modulation condition to the driven Floquet system. In the time-independent model, one requires that J j does not attain arbitrarily small values (for arbitary j), simply leading to the condition J > A J . Here, instead, we notice that the equations of motion (7) are left unchanged by the following transformation: Because of such mapping between J and J → 2π/T J − J, we define the weak modulation regime as follows It is readily seen that Eq. (8) recovers the expected condition J > A J of the undriven model when taking the appropriate limit of T → 0. On the other hand, Eq. (8) can only be satisfied for T J A J < π/4.

A. Majorana edge modes
To study the edge modes, we consider a semi-infinite chain, where it is possible to search for Γ s in the form of Majorana operators. For Floquet systems, such nontrivial edge states can lie at zero quasienergy or at the edge of the Brillouin zone ( s = π), and we denote them as Γ 0 and Γ π , respectively. To find analytic expression for the phase boundaries associated with Γ 0,π , in A we write explicitly the relevant eigenproblems in terms of the α j , β j coefficients of Eq. (6). In particular, we find that where we introduce the 2 × 1 vectors v j = (β j , α j+1 ) T and the 2 × 2 transfer matrices M j , where c j , s j and c b , s b are defined after Eq. (7). Meanwhile, the starting vector v 0 can also be evaluated explicitly as where α 0 is fixed by the normalization of Γ 0 and is irrelevant for computing the phase boundaries. When Q takes the value q n+1 /q n , the transfer matrix M j is periodic with respect to its index j with period q n , namely M j = M j+qn . The total transfer matrix for one period is therefore In B we prove that the vector v 0 is an (unnormalized) eigenvector of the transfer matrix M qn , and the corresponding eigenvalue λ 0,qn satisfies (13) In the incommensurate limit n → ∞, if λ 0 = lim n→∞ λ 0,qn is less than 1 there exists a Majorana eigenmode with zero quasienergy; otherwise, the edge mode does not exist. Therefore, the phase boundaries are found by setting the left-hand side of Eq. (13) to zero. Furthermore, in the incommensurate limit the sum in the right hand side of Eq. (13) can be approximated by an integral, giving where J(θ) = J + A J cos θ in the integrand and we set T b = T J = 1, to simplify the notation (we follow this choice from now on). In the same manner, the equation to determine the phase boundary for the π edge mode is From Eqs. (14) and (15) we calculate the phase boundaries numerically and the phase diagram is depicted in Fig. 1 for different values of A J . We observe that the parameter space is divided into distinct phases, characterized by the number of Majorana edge modes (or, equivalently, by different magnetic order of the vacuum state of fermionic excitations) In Fig. 1(a) we have A J < π/4 and the weak modulation regime [see Eq. (8)] is indicated by the hatched region. The phase boundaries determined by Eqs. (13) and (15) have a non-analytical behavior at the edges of the weak modulation region. This phase diagram coincides with that of the clean one if we take A J = 0, where the phase boundaries are two straight lines [10]. For weak modulation, all Ising couplings have the same sign and the phase transition is in the conventional Ising universality class, which is verified by the linear asymptotics of the spectrum around the transition points, see Fig. 2(a). Thus, the dynamical exponent is z = 1 in this case. On the other hand, strong modulation introduces a finite density of broken links and drives the transition unstable, making it fall into a new universality class. The numerically extracted exponent is z ≈ 1.8, see Fig. 2(a). This value is close to that estimated in Refs. [55,56] using a saturated Harris-Luck criterion for the undriven TFIM, suggesting the two are in the same QP-Ising universality class.
In Fig. 1(b) we show a representative phase diagram with A J > π/4, when the Ising coupling are always strongly modulated and the phase boundaries are more complex. Here we still find non-analytic features marking the edges of a middle region in J, where we obtain a spectrum with vanishing gaps at both 0 and π quasienergies (not shown). On the other hand, Fig. 2(a) shows that the spectrum in the weakly modulated region has gaps both around 0 and π, and only one of the two gaps can disappear at the phase boundaries or in the strongly modulated region. The latter behavior is also found in the outer region of Fig. 1(b). Based on this observation, we interpret the middle region of Fig. 1(b) as overwhelmingly modulated. For our purpose of seeking time crystals, we will mainly focus on the case A J < π/4 in the rest of the paper.

B. Localization of Excitations
The fact that the quasiperiodic Ising coupling can be approximated by a sequence of periodic potentials allows to impose an analytical upper bound for the spectral measure and to calculate numerically the sum of the bandwidth of all bulk bands in a finite chain, which we call the total bandwidth (TBW) here. The TBW essentially measures the fraction of extended states in the whole Floquet spectrum quantitatively. The analytical estimation in Ref. [56], which is also applicatable to our case, shows that in the strong modulation regime (white region for A J < π/4 and the whole parameter space for A J > π/4) the spectral measure vanishes for infinite system size, excluding the existence of extended states. In the weak modulation regime, however, we refer to numerics to identify the localization-delocalization transitions (marked by black lines in Fig. 2(b)). As expected, the states close to the lines of vanishing gap are forced to be extended, as a result of the Ising universality class, whereas in the vicinity of b = 0, π/2 all states remain localized (marked by black lines in Fig. 2(b)).
Another commonly adopted diagnostic of localization is the inverse participation ratio (IPR) calculated at different quasienergies, which gives an energy-resolved characterization of localization and is capable of detecting mobility edge. For a given Floquet eigenstate expressed in the Majorana representation, the IPR is defined by where α j (β j ) is the amplitude on even (odd) Majorana site. Finite size scaling of the IPR at different energy density, IPR ∼ 1/q α with q system size, provides an efficient method for detecting localization. For extended states, α = 1; for localized states α = 0; while for critical states this exponent lies in between 0 < α < 1. We present numerical results in Fig. 2(a) for the weak modulation regime (left two panels where J = 0.7) and the strong modulation regime (right two panels where J = 0.3) regarding the energy-resolved IPR and the exponent α of its finite size scaling relation. We see that in both cases the localization-delocalization transitions depend on energy, signaling the presence of a mobility edge. However, delocalized states in the two regimes are clearly different in that strong modulation brings all delocalized states into critical states with multifractional wavefunctions, coinciding with the theoretical prediction of vanishing spectral measure in this region.

C. Long Range Magnetic Order of Excitations
The phase diagram shown in Fig. 1(a) can also be understood in terms of magnetic order of the vacuum state, satisfying Γ s |0 = 0. By conventionally choosing s ∈ [0, π], the Γ s operators and their vacuum are uniquely defined. For simplicity, although the notion of ground state is meaningless for Floquet systems, we will sometimes refer to the eigenstates generated by Γ † s as single-particle 'excitations'. However, we emphasize that the choice of |0 has no special role in the present driven model. Instead, it is simply a reference state from which one can obtain the other Floquet eigenstates by applying the quasiparticle creation operators, and all the physical properties of the model are independent of the choice for |0 .  Fig. 1(a)). In the left two panels J = 0.7 in the weak modulation regime, while in the right two J = 0.3 in the strong modulation regime. Color in the upper panels is given by the ratio -log IPR/ log qn, a rough estimate of the finitesize scaling exponent α. Dashed lines mark how gaps close around transition points, from which we extract the dynamical exponent z. (b) TBW in the weak modulation regime for AJ = 0.6 and a chain of length qn = 144 (n = 12) with periodic boundary condition. Black solid lines in the upper panel are the numerically evaluated boundary of the two fullylocalized phases. In the lower panel, we take a line cut at J = 0.7 and show convergence of the TBW as a function of the linear size qn.
As in the case of the undriven TFIM [57], no edge states imply trivial band topology and absence of longrange order, namely a paramagnetic (PM) phase; While one edge state, no matter if it is at 0 or π quasienergy, indicates nontrivial band topology and long-range ferromagnetic (FM) order. Finally, as discussed for the clean quantum Ising chain [58] (J j = J), the new phase with two edge states, which is unique to Floquet systems, also has trivial bulk bands thus no long-range correlation in the bulk. Alternatively, in terms of spin observables, we can characterized different phases by the nature of their excited states. As we will discuss shortly below (see Fig. 3), the "0" and "π" phases eigenstates are composed of domain walls. Furthermore, due to the Z 2 symmetry (see the discussion in Sec. III D), the eigenstates form pairs with fixed quasienergy splitting (zero and π, respectively for the two phases). For this reason, we shall also call these phases FM phases. On the other hand, domain walls are absent in the other two phases so we shall call them PM phases.
A remarkable consequence of the localization of excitations in the FM phases is the preservation of longrange order in a generic "excited states", which is expressed as the Slater determinant of single-particle excitations. Nonvanishing long-range order in excited states follows from the fact that in the FM phase excitations are domain-wall-like, and a single domain wall simply revert the sign of the correlator of the vacuum state. In the localized phase, each excitation is composed of a few domain walls and is immobile. This explains what we see in Fig. 3 where the correlation function changes sign many times, corresponding to the presence of many pinned domain walls. However, the correlation function never decay to zero at long distance. Long-range order protected by localization in almost all eigenstates is crucial in view of time crystal, as it helps to build the robust oscillation of magnetization in the dynamical process.

D. Z2 Time crystal Order
The second consequence of localized excitations is the realization of time crystals. First, let us point out the major difference between the two fully localized FM phases in the vicinity of b = 0 and π/2. This is best illustrated in an open chain. Let us denote bulk Floquet excitations by Γ s , in ascending order of the quasienergy s . Then any excited Floquet eigenstate is given by |s 1 , s 2 , · · · , s j = Γ † s1 Γ † s2 · · · Γ † sj |0 hosting many excitations with |0 the aforementioned vacuum states. In the fully localized phase close to b = 0 (see Fig. 2), excited states come in nearly degenerate pairs, e.g., the pair Γ † s1 Γ † s2 · · · Γ † sj |0 and Γ † 0 Γ † s1 Γ † s2 · · · Γ † sj |0 ; On the contrary, in the fully localized phase close to b = π/2, excited states form pairs with π splitting, namely Γ † s1 Γ † s2 · · · Γ † sj |0 and Γ † π Γ † s1 Γ † s2 · · · Γ † sj |0 . This spectral pairing becomes exact in the thermodynamic limit (TDL). Remarkably, spectral pairing with a π quasienergy splitting is a prerequisite for the formation of Z 2 time crystals, where the discrete time translation symmetry is spontaneously broken and period doubling of the order parameter is observed. Combined with the long-range FM order in any Floquet excited states, we may anticipate that the fully localized FM phase close to b = π/2 is a Z 2 time crystal.
We can establish the Z 2 time crystal order by examining the dynamics of correlation functions, which are amenable to numerical calculations as they respect the Ising symmetry. The equal-time correlation ψ(nT )|σ x i σ x j |ψ(nT ) detects the presence or absence of FM order in the long run. In Fig. 4 we plot the correlator ψ(nT )|σ x i σ x j |ψ(nT ) as a function of time starting from the symmetry-unbroken ground state of a clean TFIM, for both the clean and the quasiperiodic modulated cases. For the clean model, the correlator decays exponentially to zero, indicating vanishing long-range order in the stationary state; while for the quasiperiodic model, it quickly relaxes to a nonvanishing value with some persistent fluctuation. In fact, expanding the initial state in the basis of Floquet eigenstates |ψ(0) = α,s C α,s |α, s where the index s accounts for possible degeneracy of eigenstates, we can separate the stationary value from the oscillation as [59] where σ x i σ x j D = α,s,s C * α,s C α,s α, s |σ x i σ x j |α, s denotes the generalized "diagonal ensemble" when the systematic degeneracy of Floquet eigenstates |α, s is taken into consideration, and the function f (Ω) = α =α,s ,s e i(E α ,s −Eα,s)t C * α ,s C α,s α , s |σ x i σ x j |α, s δ(Ω− E α ,s + E α,s ) is the density of states weighted by the amplitudes of the initial state. For integrable systems (like models solvable by the JW transformation), the fluctuation in the equal-time correlation function is solely determined by f (Ω) and is related to the property of the single-particle spectrum [59,60]. A continuous single-particle spectrum implies vanishing fluctuation in the long-time limit whereas a pure-point spectrum implies persistent fluctuation even in the TDL. Thus we conclude that for the quasiperiodic chain the equal-time correlator fluctuates around the diagonal ensemble value σ x i σ x j D with no decay for any system size. We have verified this statement is true (not shown here) for very long time.
The equal-time correlation function cannot distinguish the time crystal phase from a normal FM phase. This can be remedied by examining the autocorrelation function σ x j (nT )σ x j . In Fig. 4(b) and (c), we plot the autocorrelation function of the central spin for the clean and quasipeioridic chains. In both cases, we see the autocorrelation function decays at an early stage then displays revival due to finite size effect. We define the time t * at which the autocorrelation function decays to its minimum for the first time as the lifetime of the time crystal in a finite system and investigate its scaling with the system size (see insets in Fig. 4). For the clean case, the lifetime scales linearly with the size as t * ∼ q while for the quasiperiodic case it scales exponentially as ln t * ∼ q. The exponential dependence of t * on system size q supports a stable time crystal phase in the TDL.

IV. CONCLUSIONS
We have investigated the integrable TFIM subject to periodic driving and quasiperiodic modulation of the Ising coupling. Thanks to the JW transformation and the special structure of the Floquet operator, we were able to obtain analytical expressions for the Majorana edge modes. By applying the arguments in Ref. [56], we show that in the strong modulation regime the spectral measure vanishes everywhere in the parameter space, while we resort to numerics in the investigation of spectral measure in the weak modulation regime and show that two fully-localized phases exist close to the two special values b = 0, π. We analyzed two consequences of fullylocalized single-particle excitations, namely long-range FM order of eigenstates and π spectral pairing, based on which a Z 2 time crystal phase was anticipated. We also presented numerical results for relatively large systems in support of the existence of time crystals in our model. Our work will be a good starting point for future works on robustness of the time-crystalline order when integrability-breaking perturbations are included.
where we made use of the periodicity of the Ising coupling J j = J j+qn in the second line. The eigenvalue can be written explicitly as and taking the absolute value and logarithm on both sides yields Eq. (13). Note that, due to det M qn = 1 the other eigenvalue should be 1/λ 0,qn and the corresponding eigenvector can be obtained readily.