Eternal discrete time crystal beating the Heisenberg limit

A discrete time crystal (DTC) repeats itself with a rigid rhythm, mimicking a ticking clock set by the interplay between its internal structures and an external force. Discrete time crystals promise profound applications in precision timekeeping and other quantum techniques. However, it has been facing a grand challenge of thermalization. The periodic driving supplying the power may ultimately bring DTCs to thermal equilibrium and destroy their coherence. Here we show that an all-to-all interaction delivers a DTC that evades thermalization and maintains quantum coherence and quantum synchronization regardless of spatial inhomogeneities in the driving ﬁeld and the environment. Moreover, the sensitivity of this DTC scales with the total particle number to the power of 3 / 2, realizing a quantum device of measuring the driving frequency or the interaction strength beyond the Heisenberg limit. Our work paves the way for designing nonequilibrium phases with long coherence time to advance quantum metrology.


I. INTRODUCTION
A periodic driving continuously pumps energies into a discrete time crystal (DTC) [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15] and eventually heats it up to infinite temperature [16][17][18].Many schemes have been proposed to slow down thermalization [5][6][7]10], including many-body localization (MBL), Floquet prethermalization, and cryptoequilibrium.Compared with other schemes retaining the coherence of DTCs within certain timescales, MBL is of particular interest.Disorder breaks an interacting system into localized l-bits [19], suppressing thermalization up to an arbitrarily long timescale.However, most studies so far have considered homogeneous drivings.In practice, the driving field may vary across a DTC and local perturbations may further amplify spatial inhomogeneities, both preventing constituents of the DTC from synchronization and impeding applying DTCs in quantum technologies.Whereas MBL could stabilize a DTC against weak inhomogeneous perturbations to π rotations [20], it is no longer powerful in the presence of strong inhomogeneities, as exponentially decayed couplings between l-bits in MBL have readily weakened the synchronization between remote parts.Fundamental questions naturally arise: how to access a DTC that could maintain quantum coherence and synchronization in the presence of arbitrarily strong inhomogeneities of driving fields and local perturbations and how to implement such a DTC to promote the precision of quantum metrology.
In this paper we present a type of DTC that has a number of unique features distinct from previously studied ones.Discrete time crystals in the literature survive a small deviation of the driving fields from uniform π pulses.In contrast, our DTC is stable against arbitrarily strong perturbations in both homogeneous and inhomogeneous pulses.It could also start from any initial state, not necessarily a superposition of only two eigenstates of the Floquet operator.Meanwhile, the hypersensitivity of our DTC to interaction strength makes it a promising quantum device to measure interactions beyond the Heisenberg limit, unlike other DTCs not sensitive to interactions within a finite range.

II. MODEL
We consider N spins 1/2 described by a Hamiltonian H = H int + n H pul δ(t − nT ), where As shown in Fig. 1(a), J is the strength of an all-to-all interaction considered in the Lipkin-Meshkov-Glick model [21].In addition, S i = 1 2 σ i , where σ i are Pauli matrices ( h = 1), H pul represents periodic pulses applied on spins, and θ i determines the angle rotated by the ith spin about the y axis.The dependence of θ i on i characterizes the spatial inhomogeneity of pulses.In the time evolution, H int and H pul control the dynamical phases and the geometric phases, respectively, both of which depend on the numbers of spin excitations.Since H int and H pul do not commute, their interplay will lead to a perfect revival of the initial state.Equation (1) can be realized using spins 1/2 coupled to a cavity or a waveguide [22,23] or z leads to an effective π rotation about the z axis between 2nT + and (2n + 1)T − such that any initial state returns to itself after 2T for any θ.
particles with long-range interactions whose ranges are much larger than the system size.The equivalence between spins 1/2 and bosons also provides a natural realization [24].
We prove that, when JT = π , any initial state returns to itself at t = 2nT − for any N ∈ 2Z and θ i as an arbitrary function of i.Here t − (t + ) denotes the time right before (after) applying a pulse.This perfect revival delivers an eternal DTC evading thermalization and equipped with a strong synchronization even in a noisy environment.Previous works on normalized all-to-all interactions have considered the small-J limit of Eq. (1) [24], not the optimal choice of JT discussed here.
Consider an initial state with m spins up and N − m spins down | (0 − ) = i |η i , where η = ↑, ↓.After the first pulse, Here | (0 + ) becomes a superposition of 2 N states.Each state is obtained from flipping s spins up and k spins down of | (0 − ) , as shown in Fig. 1(b).Each state acquires a dynamical phase e −iϕ 1 imposed by H int from t = 0 + to t = T − .The second pulse flips spins again, followed by H int imposing another dynamical phase e −iϕ 2 from t = T + to t = 2T − .We define To return to | (0 − ) , the s (k) spins up (down) flipped by the first pulse need to flip back to spins up (down) during the second pulse; 2 N such pathways allow the system to return to | (0 − ) .The contribution to A from each pathway is (−1) k+s j∈ F cos 2 ( , where (−1) k+s comes from flipping k + s spins 1/2 twice, equivalent to the geometric phase from rotating these spins about the y axis for 2π .Here F ( F ) denotes the collection of flipped (unflipped) spins.As each of these 2 N states is an eigenstate of H int , ϕ 1 = (m − s + k)(m − s + k − N )π and ϕ 2 = m(m − N )π when JT = π .The m-independent terms have been dropped.Since N ∈ 2Z and e iZ 2 π = (−1) Z ∀ Z ∈ Z, we have e −i(ϕ 1 +ϕ 2 ) = e iπ{2[m 2 +m(k−s−N )−ks]+N (s−k)+k 2 +s 2 } = (−1) k+s , which cancels the previously obtained geometric phase, and thus A = 2 ).Since F denotes the sum over all 2 N choices of flipping the N spins in | (0 − ) , This discussion applies to any initial product state and any t ∈ [2nT − , 2(n + 1)T − ].Thus, any initial state returns to itself at t = 2nT − .Unlike traditional spin-echo schemes using tailored pulses [25], we implement interactions, one source of the decoherence, to overcome the other, the inhomogeneities, so as to access a perfect dynamical localization, an analogy to the Anderson localization in the Hilbert space [26].This interaction-induced spin echo applies to a broad class of systems to extend the coherence time.
For spatially uniform pulses, a simpler proof exists.Here H is rewritten as where L = i S i .Equation ( 5) is equivalent to the kicked top model for spin L [27], where L = N 2 .The propagator from t = 2nT − to t = 2(n + 1)T − is U JT (2T ) = e −iJT L 2 z e −iθL y e −iJT L 2 z e −iθL y .As e −iπL 2 z = e −iπL z applies to any integer L (or even N), U π (2T ) = e −2iπL z e e iπLz (−iθL y )e −iπLz e −iθL y = 1.As shown in Fig. 1(c), any state on the Bloch sphere returns to the original place after 2T .
z and e −iπL z are no longer identical, and such a DTC with a period of 2T does not exist.In contrast, if we consider spin 1 in Eq. ( 1), such an even-odd effect is absent, since L ∈ Z ∀ N.
Here U π (2T ) = 1 means that the quasienergy spectrum of H eff , where U JT (2T ) = e −iH eff , has 2 N degenerate eigenstates.Whereas this looks similar to the noninteracting case when θ i = π , a conceptual difference is that the degeneracy here is stable against any perturbations in θ i , unlike noninteracting systems, where any infinitesimal deviation from a homogeneous π pulse lifts the degeneracy, breaks the integrability, and suppresses DTCs.Similar to other models, the period doubling comes from the spontaneous time-translation symmetry breaking [5,6,20,24].When JT = π , the Floquet eigen- Comparison between the all-to-all interaction and a power-law potential with α = 3.Here N = 14.(a)-(d) Uniform pulses, with w s = 0.The DTC with all-to-all interactions (dots) is unaffected by the deviation .With the power-law potential (curves), increasing leads to the suppression of P(2nT ) and M z (2nT ) and the growth of S(2nT ) and Q(2nT ).(e)-(h) Keeping θ = π and increasing the spatial inhomogeneities w s , the DTC with the power-law potential is suppressed.The DTC with all-to-all interactions remains stable.See Appendix B for more results of θ = π .
Rabi oscillation between the two Floquet eigenstates leads to a period of 2T .Previous works have mainly focused on uniform π pulses, where the Floquet eigenstate 1 is a special case of our results.In particular, results here apply to any uniform and nonuniform pulses and thus are far more general.

III. STABILITY AGAINST SPATIAL INHOMOGENEITIES
We compare our model to the power-law interaction H = H int + n H pul δ(t − nT ), where Starting from | (0 − ) = i |↑ i , we compute Here P(2nT − ) characterizes the quantum memory of the initial state, M z (2nT − ) denotes the z component of the total spin, captures the absorption of energy, and S(2nT − ) is the halfchain entanglement entropy.
When θ i = θ , a finite J in Eq. ( 6) restores the quantum coherence, if = θ − π is small [7][8][9].With increasing , characterizes the absorption of the energy.In addition, E ∞ = 2 −N j j|H int | j is the energy at infinite temperature and j sums over all eigenstates of H int .These results signify the thermalization at large .Taking into account spatial inhomogeneities, as shown in Figs.2(e)-2(h), we choose a random θ i from [ θ − w s , θ + w s ] with a constant probability.When w s is finite, the thermalization becomes even faster and Q approaches 1, indicating thermalization to infinite temperature.For power-law interactions, dynamical phases controlled by interactions cannot cancel geometric phases induced by pulses.It is impossible to obtain a constructive interference between all pathways.Whereas on-site disorders are introduced to create MBL to slow down the thermalization [7][8][9], an intrinsic drawback is that MBL weakens the synchronization when the pulses are nonuniform (Appendix A).In contrast, P(2nT − ) and M z (2nT − ) of the all-to-all interaction are unaffected by w s , and Q(2nT − ) and S(2nT − ) remain zero, directly reflecting the robustness of this eternal DTC.

IV. APPLICATIONS
The perfect revival at t = 2nT − comes from the same phase of all 2 N pathways of returning to | (0) when JT = π .Once JT = π , the larger N is, the more rapidly the phase varies with changing pathways.In the large-N limit, this DTC becomes supersensitive to JT and serves as a high-precision device to measure either J or T .
It is time consuming to solve more than 14 lattice sites using exact diagonalization when inhomogeneities exist.We focus on homogeneous systems.It is expected that the lower bound of the results of an inhomogeneous distribution θ i ∈ [ θ − w s , θ + w s ] could be estimated using homogeneous θ i = θ ± w s .As an example, we consider θ i fixed at π/4.As shown in Fig. 3(a), P(2nT − ) quickly vanishes if |δ| π/N 3/2 , where δ = JT − π .Whereas using highly entangled states in linear interferometry could beat the standard quantum limit 1/ √ N and access the Heisenberg limit 1/N, nonlinear interferometry using k-body interactions could beat the Heisenberg limit and reach the ultimate limit N −k [28][29][30][31][32][33].Our DTC represents a different category of nonlinear metrology using periodic drivings to beat the Heisenberg limit.
The P(2T ) captures short-time dynamics.Figure 3(c) shows that the dependence of P(2T ) on JT has a narrow peak centered at π , whose width ∼1/N 3/2 .Such scaling can be obtained analytically (Appendix C 1) and is verified numerically [inset of Fig. 3(c  response at half of the frequency of the periodic driving.Its dependence on JT also has a peak around π .We define the full width at half maximum as JT and find both numerically and analytically that JT is proportional to 1/N 3/2 (Appendix C 1).
The DTCs discussed in the literature are stable within a finite range of both the interaction strength and a uniform deviation of θ i from π .In contrast, our DTC is stable against any spatial fluctuations in θ i and while being supersensitive to JT .In practice, it is easier to control J and T than the N local parameters θ i in a noisy environment, where θ i may not have any correlations at different locations.Moreover, our DTC could measure JT with high precision beyond the Heisenberg limit.It mimics a supersensitive clock.If the frequency of the external field ω d = 2π/T is fixed, J corresponds to some internal parameter of a clock (for instance, the length of a pendulum clock) and needs to be tuned with a precision of 1/N 3/2 to deliver rigid ticks at t = 2nT .Otherwise, as shown in Figs.3(c) and 3(e), once JT is beyond the length 1/N 3/2 centered at JT = π , both P(2T ) and P( f ) quickly decrease and the DTC stalls to avoid errors in the timekeeping.From JT = π , the precision of J can be estimated as , where d /ω d characterizes the precision of the driving frequency.When , the uncertainty of J approaches the precision limit of d /ω d .Whereas the precision of d /ω d is up to 10 −19 in the terahertz regime [35], typical experiments on ultracold atoms, ion traps, and nitrogenvacancy centers have interaction strengths ∼10 2 -10 5 Hz.In such a regime, the precision of d /ω d could be 10 −6 and higher.Our results thus provide an application of precision timekeeping in many-body physics.
Fixing J, our DTC could also gauge the frequency.Only when T deviates from π/J within 1/N 3/2 could a driving field induce a long-lasting dynamics.Unlike atomic clocks using a narrow linewidth transition, here many-body effects determine the driving frequency.The rotated angle θ is arbitrary such that this DTC functions in a nonideal environment, unlike previous works requiring precise control of pulses in nonlinear metrology without periodic driving [29,30,36].
Previous conclusions apply to a generic long-range interaction if its range is much larger than the size of the system.For instance, with decreasing α, the range of a power-law potential in Eq. ( 6) increases.When α = 0, it is equivalent to the all-to-all interaction.Figure 4 shows results for N = 14.With decreasing α down to zero, P(2nT ) and M z (2nT ) increase and approach the result of the all-to-all interaction.A small α = 0.04 readily provides us with a good approximation of the all-to-all interaction in such a finite system.
Both interactions and external drivings are crucial for DTCs.We hope that our work will stimulate more studies and N 2 = L is the total angular momentum As discussed in the main text, e −iπL 2 z = e −iπL z is satisfied for any even particle number N. When δ is small, e −iδL 2 z can be written as We thus obtain Noting that we obtain where α ≡ tan 2 θ 2 and z by e −iπ L 2 z , we have ignored the expansion of the wave function in the latitude direction that gives rise to a high-order correction to P(2nT ) at small times.is a coherent state pointing along θ and φ.In the large-N limit, which represents a narrow Gaussian centered at k = 0. Substituting I (θ, k) in Eq. (C5) using Eq.(C11), we obtain Eq. (C1).

Scalings of P( 1 2T )
We have also obtained an analytical form for P( f ), the Fourier transform of P(2nT ).As shown in Fig. 8, starting from an initial state at the north pole, the state at t = 2T − covers a finite small region near the north pole if δ = JT − π is small.The length scales of the longitude and latitude directions are proportional to δ and δ 2 , the latter of which can be ignored in the small-δ limit.Thus, we make use of the following approximation to capture the dynamics in the small-δ limit: Here P(2nT ) is written as Using the identities e −iπL z e −iθL y e −iπL z = e iθL y and e −iπL 2 z = e −iπL z , Eq. (C17) can be written as Applying Eq. (C1), we obtain Equation (C19) recovers Eq. (C1) when n = 1.As shown in Fig. 7(a), this expression captures well the initial decay of P(2nT ).However, it cannot describe the revival of P(2nT ) at later times for certain JT due to the approximation made in Eq. (C16).The power spectrum is therefore written as where M is the cutoff required in numerics.In the last step, we have used the fact that, for small n, |ψ (2nT + T ) is located at a place on the Bloch sphere away from the north pole, provided θ is not small, and thus P(2nT When nNδ 1 and θ = 0, π/2, π, Eq. (C19) becomes P(2nT ) = e −n 2 sin 2 (2θ )δ 2 N 3 /16 and Eq.(C20) is rewritten as P(1/2T ) ≈ 2 √ π erf 1 8 δMN 3/2 sin 2θ sin(2θ )δMN 3/2  .(C21) In the limit M → ∞, P(1/2T ) → 1
The comparison between this analytical result and the numerical one is shown in Fig. 3(e).
Here M z (2nT ) and Mz (1/2T ) do not have simple analytical forms.We have numerically evaluated them and the scaling of Mz (1/2T ) with N is shown in Figs.

APPENDIX D: QUANTUM FISHER INFORMATION
When JT = π , the quantum Fisher information is written as )

APPENDIX E: NUMERICAL SIMULATION METHODS
In this Appendix we present some details of the numerical simulations used to produce Figs.2, 3, and 7. We write the time-dependent many-body wave function | (t ) as a superposition of Fock states, which are eigenstates of S z .In the presence of H int , every Fock state acquires a dynamical phase factor.Then H pul is applied to flip the spins at t = nT to obtain the evolution of the many-body wave function.
The entanglement entropy S in Fig. 2 is calculated by first tracing out half of the system and obtaining a density matrix ρ B for the other half of the system.We then diagonalize ρ B and obtain all its eigenvalues v i .Then S is calculated by using S = − i v i ln(v i ).
The power spectra shown in Figs. 3 and 7 are obtained as follows.First P(t ) and M z (t ) are evaluated in the time interval between t = 0 and t = MT , where M represents the longest time we consider in the simulation.We then use P( f ) = 1 M M−1 n=0 e i2πnT f P(nT ) to calculate the power spectra.We use the full width at half maximum to characterize the width of power spectrum.Namely, when the value of P or M z is decreased to half of its maximum, we define twice the deviation of JT from π as the width.

APPENDIX F: EXPERIMENTAL REALIZATION
In this Appendix we discuss the experimental realization of our model.Whereas we used idealized δ-like kicks in the discussions in the main text, our results can be straightforwardly generalized to pulses with finite widths.If we consider the Hamiltonian where τ represents the finite width, we only need to replace JT = π by J (T − τ ) = π and other results remain unchanged.In a finite system, as shown in the main text, to qualitatively demonstrate our results for a system consisting of 14 spins, a power-law interaction width α = 0.04 readily provides us with a good approximation of the all-to-all interaction.In fact, the only requirement is that the range of interaction is much larger than the size of the system.Meanwhile, the all-to-all interaction can be accessed by using photons in cavities or waveguides to couple atoms at different locations.
We consider the Hamiltonian where the first term represents a single-photon mode and the second term denotes the interaction between photons and atoms.We have assumed that the local Hamiltonian acting on each spin vanishes.In addition, S ± i = S x i ± S y i .Without the second term, the eigenstates are |n i |η i , representing n photons and a Fock state of atoms, and η = ↑, ↓.Consider two specific atoms at sites i = j; their couplings induced by the photon can be derived using second-order perturbation.For instance, the off-diagonal term is written as This is equivalent to adding a term − 2 hω S + j S − i to the unperturbed Hamiltonian.The diagonal couplings can be obtained by a similar means.The full effective Hamiltonian after eliminating the photon mode is written as

(F4)
Since i S i is conserved, the first two terms commute with the last term and the last term exactly matches an all-to-all interaction; we thus created an all-to-all interaction with an equivalent strength J = 2 2 hω .Whereas the above discussion does not require the leakage of photons from the cavity, a "bad" cavity with leaking photons has the unique advantage of suppressing the heating that may be caused by the driving [37].Thus, our results can be generalized to the full model including both the atoms and photons.
Whereas the above scheme is relevant to small systems, in which a fine-tuning of the local effective magnetic field is doable, our model can also be implemented by a two-mode bosonic system, in which the all-to-all interaction naturally exists.For instance, we consider a Bose-Einstein condensate in a double-well potential H = g 1 n l (n l − 1) + n r (n r − 1) 2 + g 2 n l n r + n (θ a † l a r + H.c.)δ(t − nT ), (F5) where g 1 and g 2 represent on-site and intersite interactions, respectively, and θ denotes the tunneling.If we map the left and the right site to the spin up and the spin down, respectively, J = g 1 + g 2 directly corresponds to the all-to-all interaction, and the tunneling term is mapped to H pul .Such a mapping allows us to implement all results in the spin model to a large bosonic system with 10 4 particles and more.

FIG. 1 .
FIG. 1.(a) All-to-all interactions (brown arrows) between spin 1/2 (blue spheres attached to arrows).(b) Perfect revival of an arbitrary initial state due to constructive interference among all pathways.Dashed and solid boxes highlight the k spins up and s spins down flipped by the first pulse, leading to geometric phases (−1) k and (−1) s , respectively.Triangles on the time axis represent H pul .(c) Rotations of a spin L (yellow arrow) on the Bloch sphere.When JT = π , the nonlinear term JL 2z leads to an effective π rotation about the z axis between 2nT + and (2n + 1)T − such that any initial state returns to itself after 2T for any θ.

FIG. 3 .
FIG. 3. Sensitivity to JT .Curves (dots) are numerical (analytical) results.Orange, blue, and green colors in (c)-(f) represent N = 100, 200, and 400, respectively.(a) P(2nT ) and (b) M z (2nT ) as functions of n at various JT .When |JT − π | π/N 3/2 , both quantities quickly decrease down to zero.Here N = 200 has been used in the calculation.(c) P(2T ) and (d) M z (2T ) as a function of JT .For a fixed N, P(2T ) and M z (2T ) have narrow peaks centered at JT = π .In contrast, M z (2T ) has an additional fast oscillation; the dashed curve highlights the analytical result of its profile, whose width is denoted by black arrows.(e) P(1/2T ) and (f) Mz (1/2T ) are also featured by narrow peaks around JT = π .Here M = 200 is used in numerics.Insets in (c)-(f) show scalings of widths of peaks with N. (g) Quantum Fisher information I π (2nT ) as a function of n.(h) I π (2nT ) is proportional to N 3 .The value θ i = π/4 is used in all panels.

FIG. 4 .FIG. 6 .
FIG. 4. Power-law potentials for N = 14.(a) and (b) Uniform rotations with θ = 0.95π .Here α = 0 corresponds to the all-to-all interaction.With decreasing α, the results of the power-law potentials approach those of the all-to-all interaction.Other parameters are the same in (a) and (b).(c) and (d) Inhomogeneous rotations with θ = π and w s = 0.1π .

FIG. 7 .
FIG. 7. Sensitivity of the DTC to JT when θ = π/2.Dots are analytical results and curves are the numerical results.Orange, blue, and green colors in (c)-(f) represent N = 100, 200, and 400, respectively.(a) P(2nT ) and (b) M z (2nT ) as functions of n for various JT .When |JT − π | π/N (here N = 200), both quantities quickly decrease down to zero.(c) P(2T ) and (d) M z (2T ) at a fixed time t = 2T as functions of JT .For a fixed N, both quantities are featured with narrow peaks centered at JT = π .Insets show the scalings of the widths (full width at half maximum) of the peaks with N. Power spectra (e) P(1/2T ) and (f) Mz (1/2T ) are also featured with narrow peaks around JT = π .Whereas they exhibit nonmonotonic behaviors near JT = π , both quantities vanish when |JT − π | π/N.Insets show the scalings of the widths of the peaks with N. (g) Quantum Fisher information I π (2nT ) as a function of n.(h) I π (2nT ) is proportional to N 2 .The value θ i = π/2 is used in all panels.