Time crystals in the driven transverse field Ising model under quasiperiodic modulation

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

$$\begin{equation}{U}_{\text{F}}=\mathrm{exp}\left[-\text{i}{T}_{J}\sum _{j}{J}_{j}{\sigma }_{j}^{x}{\sigma }_{j+1}^{x}\right]\mathrm{exp}\left[-\text{i}{T}_{b}b\sum _{j}{\sigma }_{j}^{z}\right].\end{equation} \tag{ 3 }$$

There are several symmetries of this Floquet operator. The Z₂ Ising symmetry, denoted by the symmetry operator $P={\prod }_{j\enspace }{\sigma }_{j}^{z}$, is inherited from the Hamiltonian (1) and satisfies PU_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 ${\sigma }_{j}^{z}\to -{\sigma }_{j}^{z}$. Instead, the T_J J reflection can be combined with φ_J → φ_J + π and local rotations on the even (or odd) sites, such that ${\sigma }_{j}^{x}\to -{\sigma }_{j}^{x}$. 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 + π.

To describe the time evolution in terms of noninteracting fermions, we introduce a pair of Majorana fermions for each spin-1/2

$$\begin{equation}{\gamma }_{2j}=\left(\prod _{i{< }j}{\sigma }_{i}^{z}\right){\sigma }_{j}^{x},\hspace{25.0pt}{\gamma }_{2j+1}=\left(\prod _{i{< }j}{\sigma }_{i}^{z}\right){\sigma }_{j}^{y}.\end{equation} \tag{ 4 }$$

The Majorana operators satisfy the anticommutation relation {γ_i , γ_j } = 2δ_ij . Then, the Floquet operator can be expressed as U_F = U₁ U₂, where:

$$\begin{align}\hfill {U}_{1}& =\mathrm{exp}\left(-{T}_{J}\sum _{j}{J}_{j}{\gamma }_{2j+1}{\gamma }_{2j+2}\right)\hfill \\ \hfill {U}_{2}& =\mathrm{exp}\left(-{T}_{b}b\sum _{j}{\gamma }_{2j}{\gamma }_{2j+1}\right).\hfill \end{align} \tag{ 5 }$$

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:

$$\begin{equation}{{\Gamma}}_{s}=\sum _{j=0}^{\infty }{\alpha }_{j}{\gamma }_{2j}+\sum _{j=0}^{\infty }{\beta }_{j}{\gamma }_{2j+1},\end{equation} \tag{ 6 }$$

and satisfy ${U}_{\text{F}}^{{\dagger}}{{\Gamma}}_{s}{U}_{\text{F}}={\text{e}}^{-\text{i}{{\epsilon}}_{s}}{{\Gamma}}_{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 (${{\Gamma}}_{s}^{{\dagger}}$) are fermionic annihilation (creation) operators, satisfying $\left\{{{\Gamma}}_{s},{{\Gamma}}_{{s}^{\prime }}^{{\dagger}}\right\}={\delta }_{s,{s}^{\prime }}$. Since Γ_s , ${{\Gamma}}_{{s}^{\prime }}^{{\dagger}}$ 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 equation (A.1)) and, with j in the bulk, yields:

$$\begin{equation}\begin{aligned}\hfill {\alpha }_{j}{c}_{j-1}{c}_{b}+{\alpha }_{j+1}{s}_{j}{s}_{b}+{\beta }_{j}{c}_{j}{s}_{b}-{\beta }_{j-1}{s}_{j-1}{c}_{b}={\text{e}}^{-\text{i}{{\epsilon}}_{s}}{\alpha }_{j},\\ \hfill -{\alpha }_{j}{c}_{j-1}{s}_{b}+{\alpha }_{j+1}{s}_{j}{c}_{b}+{\beta }_{j}{c}_{j}{c}_{b}+{\beta }_{j-1}{s}_{j-1}{s}_{b}={\text{e}}^{-\text{i}{{\epsilon}}_{s}}{\beta }_{j},\end{aligned}\end{equation} \tag{ 7 }$$

where we defined c_j =  cos  2T_J J_j , s_j =  sin  2T_J J_j and c_b =  cos  2T_b b, s_b =  sin  2T_b b.

Finally, based on equation (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: J → 2π/T_J − J, ϕ_J → ϕ_j + π, → + π, while α_j → (−1)^j α_j and β_j → (−1)^j β_j . Because of such mapping between J and J → 2π/T_J − J, we define the weak modulation regime as follows:

$$\begin{equation}{A}_{J}{< }J{< }\frac{2\pi }{{T}_{J}}-{A}_{J}.\end{equation} \tag{ 8 }$$

It is readily seen that equation (8) recovers the expected condition J > A_J of the undriven model when taking the appropriate limit of T → 0. On the other hand, equation (8) can only be satisfied for T_J A_J < π/4.

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 Γ₀ and Γ_π , respectively. To find analytic expression for the phase boundaries associated with Γ_0,π , in appendix A we write explicitly the relevant eigenproblems in terms of the α_j , β_j coefficients of equation (6). In particular, we find that ${U}_{\text{F}}^{{\dagger}}{{\Gamma}}_{0}{U}_{\text{F}}={{\Gamma}}_{0}$ is equivalent to the recurrence relation

$$\begin{equation}{\mathbf{v}}_{j}={M}_{j}{\mathbf{v}}_{j-1},\end{equation} \tag{ 9 }$$

where we introduce the 2 × 1 vectors ${\mathbf{v}}_{j}={\left({\beta }_{j},{\alpha }_{j+1}\right)}^{\text{T}}$ and the 2 × 2 transfer matrices M_j ,

$$\begin{equation}{M}_{j}=\frac{1}{{s}_{b}}\left[\begin{matrix}\hfill {s}_{j-1}\hfill & \hfill {c}_{b}-{c}_{j-1}\hfill \\ \hfill \frac{{s}_{j-1}\left({c}_{b}-{c}_{j}\right)}{{s}_{j}}\hfill & \hfill \frac{1+{c}_{j-1}{c}_{j}-{c}_{b}\left({c}_{j-1}+{c}_{j}\right)}{{s}_{j}}\hfill \end{matrix}\right],\end{equation} \tag{ 10 }$$

where c_j , s_j and c_b , s_b are defined after equation (7). Meanwhile, the starting vector v₀ can also be evaluated explicitly as

$$\begin{equation}{\mathbf{v}}_{0}={\alpha }_{0}\frac{1-\mathrm{cos}\enspace 2{T}_{b}b}{\mathrm{sin}\enspace 2{T}_{b}b}\left[\begin{matrix}\hfill -1\hfill \\ \hfill \frac{1+\mathrm{cos}\enspace 2{T}_{J}{J}_{0}}{\mathrm{sin}\enspace 2{T}_{J}{J}_{0}}\hfill \end{matrix}\right],\end{equation} \tag{ 11 }$$

where α₀ is fixed by the normalization of Γ₀ 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+{q}_{n}}$. The total transfer matrix for one period is therefore

$$\begin{equation}{\mathbf{M}}_{{q}_{n}}=\prod _{j=1}^{{q}_{n}}{M}_{j}.\end{equation} \tag{ 12 }$$

In appendix B we prove that the vector v₀ is an (unnormalized) eigenvector of the transfer matrix ${\mathbf{M}}_{{q}_{n}}$, and the corresponding eigenvalue ${\lambda }_{0,{q}_{n}}$ satisfies:

$$\begin{equation}\mathrm{ln}\vert {\lambda }_{0,{q}_{n}}\vert ={q}_{n}\left(\mathrm{ln}\left\vert \mathrm{tan}\enspace {T}_{b}b\right\vert +\frac{1}{{q}_{n}}\sum _{j=0}^{{q}_{n}-1}\mathrm{ln}\left\vert \mathrm{cot}\enspace {T}_{J}{J}_{j}\right\vert \right).\end{equation} \tag{ 13 }$$

In the incommensurate limit n → ∞, if ${\lambda }_{0}=\underset{n\to \infty }{\mathrm{lim}}\enspace {\lambda }_{0,{q}_{n}}$ 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 equation (13) to zero. Furthermore, in the incommensurate limit the sum in the right-hand side of equation (13) can be approximated by an integral, giving:

$$\begin{equation}-\mathrm{ln}\left\vert \mathrm{tan}\enspace b\right\vert ={\int }_{0}^{2\pi }\frac{\mathrm{d}\theta }{2\pi }\enspace \mathrm{ln}\left\vert \mathrm{cot}\enspace J\left(\theta \right)\right\vert ,\end{equation} \tag{ 14 }$$

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

$$\begin{equation}-\mathrm{ln}\left\vert \mathrm{cot}\enspace b\right\vert ={\int }_{0}^{2\pi }\frac{\mathrm{d}\theta }{2\pi }\enspace \mathrm{ln}\left\vert \mathrm{cot}\enspace J\left(\theta \right)\right\vert .\end{equation} \tag{ 15 }$$

From equations (14) and (15) we calculate the phase boundaries numerically and the phase diagram is depicted in figure 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)

Zoom In Zoom Out Reset image size

In figure 1(a) we have A_J < π/4 and the weak modulation regime (see equation (8)) is indicated by the hatched region. The phase boundaries determined by equations (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 figure 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 figure 2(a). This value is close to that estimated in references [56, 57] using a saturated Harris-Luck criterion for the undriven TFIM, suggesting the two are in the same QP-Ising universality class.

Zoom In Zoom Out Reset image size

In figure 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, figure 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 figure 1(b). Based on this observation, we interpret the middle region of figure 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.

The fact that the QP 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 reference [57], 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 figure 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 figure 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

$$\begin{equation}\text{IPR}=\sum _{j}\left(\vert {\alpha }_{j}{\vert }^{4}+\vert {\beta }_{j}{\vert }^{4}\right),\end{equation} \tag{ 16 }$$

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 figure 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.

The phase diagram shown in figure 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 ${{\Gamma}}_{s}^{{\dagger}}$ 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⟩.

As in the case of the undriven TFIM [58], no edge states imply trivial band topology and absence of long-range 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 [59] (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 figure 3), the '0' and 'π' phases eigenstates are composed of domain walls. Furthermore, due to the Z₂ symmetry (see the discussion in section 3.4), 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.

Zoom In Zoom Out Reset image size

A remarkable consequence of the localization of excitations in the FM phases is the preservation of long-range 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 figure 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.

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 $\vert {s}_{1},{s}_{2},\dots ,{s}_{j}\rangle ={{\Gamma}}_{{s}_{1}}^{{\dagger}}{{\Gamma}}_{{s}_{2}}^{{\dagger}}\cdots {{\Gamma}}_{{s}_{j}}^{{\dagger}}\vert 0\rangle$ hosting many excitations with |0⟩ the aforementioned vacuum states. In the fully localized phase close to b = 0 (see figure 2), excited states come in nearly degenerate pairs, e.g. the pair ${{\Gamma}}_{{s}_{1}}^{{\dagger}}{{\Gamma}}_{{s}_{2}}^{{\dagger}}\cdots {{\Gamma}}_{{s}_{j}}^{{\dagger}}\vert 0\rangle$ and ${{\Gamma}}_{0}^{{\dagger}}{{\Gamma}}_{{s}_{1}}^{{\dagger}}{{\Gamma}}_{{s}_{2}}^{{\dagger}}\cdots {{\Gamma}}_{{s}_{j}}^{{\dagger}}\vert 0\rangle$; On the contrary, in the fully localized phase close to b = π/2, excited states form pairs with π splitting, namely ${{\Gamma}}_{{s}_{1}}^{{\dagger}}{{\Gamma}}_{{s}_{2}}^{{\dagger}}\cdots {{\Gamma}}_{{s}_{j}}^{{\dagger}}\vert 0\rangle$ and ${{\Gamma}}_{\pi }^{{\dagger}}{{\Gamma}}_{{s}_{1}}^{{\dagger}}{{\Gamma}}_{{s}_{2}}^{{\dagger}}\cdots {{\Gamma}}_{{s}_{j}}^{{\dagger}}\vert 0\rangle$. This spectral pairing becomes exact in the thermaldynamic limit (TDL). Remarkably, spectral pairing with a π quasienergy splitting is a prerequisite for the formation of Z₂ 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₂ time crystal.

We can establish the Z₂ 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 $\langle \psi \left(nT\right)\vert {\sigma }_{i}^{x}{\sigma }_{j}^{x}\vert \psi \left(nT\right)\rangle$ detects the presence or absence of FM order in the long run. In figure 4 we plot the correlator $\langle \psi \left(nT\right)\vert {\sigma }_{i}^{x}{\sigma }_{j}^{x}\vert \psi \left(nT\right)\rangle$ as a function of time starting from the symmetry-unbroken ground state of a clean TFIM, for both the clean and the QP modulated cases. For the clean model, the correlator decays exponentially to zero, indicating vanishing long-range order in the stationary state; while for the QP 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 [60]

$$\begin{equation}\langle {\sigma }_{i}^{x}{\sigma }_{j}^{x}\rangle \left(nT\right)={\langle {\sigma }_{i}^{x}{\sigma }_{j}^{x}\rangle }_{D}+\int \mathrm{d}{\Omega}\enspace {\text{e}}^{-\text{i}n{\Omega}T}f\left({\Omega}\right),\end{equation} \tag{ 17 }$$

where ${\langle {\sigma }_{i}^{x}{\sigma }_{j}^{x}\rangle }_{D}={\sum }_{\alpha ,s,{s}^{\prime }}\enspace {C}_{\alpha ,{s}^{\prime }}^{{\ast}}{C}_{\alpha ,s}\langle \alpha ,{s}^{\prime }\vert {\sigma }_{i}^{x}{\sigma }_{j}^{x}\vert \alpha ,s\rangle$ denotes the generalized 'diagonal ensemble' when the systematic degeneracy of Floquet eigenstates |α, s⟩ is taken into consideration, and the function $f \left({\Omega}\right)={\sum }_{{\alpha }^{\prime }\ne \alpha ,{s}^{\prime },s}\enspace {\text{e}}^{\text{i}\left({E}_{{\alpha }^{\prime },{s}^{\prime }}-{E}_{\alpha ,s}\right)t}{C}_{{\alpha }^{\prime },{s}^{\prime }}^{{\ast}}{C}_{\alpha ,s}\langle {\alpha }^{\prime },{s}^{\prime }\vert {\sigma }_{i}^{x}{\sigma }_{j}^{x}\vert \alpha ,s\rangle \delta \left({\Omega}-{E}_{{\alpha }^{\prime },{s}^{\prime }}+{E}_{\alpha ,s}\right)$ 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 [60, 61]. 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 QP chain the equal-time correlator fluctuates around the diagonal ensemble value ${\langle {\sigma }_{i}^{x}{\sigma }_{j}^{x}\rangle }_{D}$ with no decay for any system size. We have verified this statement is true (not shown here) for very long time.

Zoom In Zoom Out Reset image size

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 $\langle {\sigma }_{j}^{x}\left(nT\right){\sigma }_{j}^{x}\rangle$. In figures 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 figure 4). For the clean case, the lifetime scales linearly with the size as t_* ∼ q while for the QP 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.