The isolated Heisenberg magnet as a quantum time crystal

TL;DR

This paper shows that the isolated Heisenberg XXZ spin chain can act as a genuine quantum time crystal, exhibiting persistent oscillations due to dynamical symmetries, without external driving or environment coupling.

Contribution

It analytically and numerically demonstrates the existence of a quantum time crystal in an isolated Heisenberg magnet, highlighting the role of dynamical symmetries and their fractal dependence on anisotropy.

Findings

The Heisenberg XXZ chain does not relax to stationarity.

Persistent oscillations are caused by dynamical symmetries.

Frequency of oscillations varies fractally with anisotropy.

Abstract

We demonstrate analytically and numerically that the paradigmatic model of quantum magnetism, the Heisenberg XXZ spin chain, does not relax to stationarity and hence constitutes a genuine time crystal that does not rely on external driving or coupling to an environment. We trace this phenomenon to the existence of extensive dynamical symmetries and find their frequency to be a no-where continuous (fractal) function of the anisotropy parameter of the chain. We discuss how the ensuing persistent oscillations that violate one of the most fundamental laws of physics could be observed experimentally and identify potential metrological applications.

Equations

$$[H,Y]=\omega Y,$$

$$\begin{array}{rl}& H=J[(\sum _j{s}_j^x{s}_{j+1}^x+s_j^y{s}_{j+1}^y+\Delta s_j^z{s}_{j+1}^z+\\ & +\alpha (\sum _j{s}_j^x{s}_{j+2}^x+s_j^y{s}_{j+2}^y+\Delta s_j^z{s}_{j+2}^z)]+\sum _{j}h{s}_j^z,\end{array}$$

$${\mathcal{M}}_{\omega }(Y_l)=\exp (\mathrm{i}2\pi H{\omega }^{-1})Y_l\exp (-\mathrm{i}2\pi H{\omega }^{-1}),$$

$${\rho }_{tGGE}=\exp (-\sum _j{\mu }_j{X}_j-{\mu }_Y(t)Y-{\bar{\mu }}_Y(t)Y^{†}).$$

$$\begin{array}{rl}⟨\psi |Y_k(t)|\psi ⟩& =\mathrm{tr}(Y_k(t){\rho }_{tGGE}),\\ ⟨\psi |X_k|\psi ⟩& =\mathrm{tr}(X_k{\rho }_{tGGE}),\end{array}$$

$$\mathrm{tr}(Y_k^l{Y}_{k^{\prime }}^{l^{\prime }†})$$

$$\mathrm{tr}(Y_k^l{Y}_{k^{\prime }}^{l^{\prime }})$$

$$O(t)\underset{t\to \infty }{=}{\alpha }_Y^O\exp (\mathrm{i}\omega t)Y+h.c.+\sum _j{\alpha }_j^O{X}_j,$$

$$O_k(t)=\int \mathrm{d}\lambda (\exp (\mathrm{i}\omega t)f(\lambda )Y(\lambda )+\exp (-\mathrm{i}\omega t)\bar{f}(\lambda )Y^{†}(\lambda )).$$

$$⟨O_k{Y}^{†}({\lambda }^{\prime })⟩=\int \mathrm{d}\lambda f(\lambda )⟨Y(\lambda )Y^{†}({\lambda }^{\prime })⟩$$

$$|\psi ⟩=2^{-N/2}(|\uparrow ⟩+|\downarrow ⟩{)}^{\otimes n},$$

$$L(\lambda )=\frac{1}{\sinh \eta }\left[\begin{array}{cc}\sin (\lambda +\eta {\mathbf{S}}^z)& \sin (\eta ){\mathbf{S}}^{-}\\ \sin (\eta ){\mathbf{S}}^{+}& \sin (\lambda -\eta {\mathbf{S}}^z)\end{array}\right].$$

$$[{\mathbf{S}}^{+},{\mathbf{S}}^{-}]$$

$$q^{2{\mathbf{S}}^z}{\mathbf{S}}^{\pm }$$

$${\mathbf{S}}^z=\sum _{n=0}^{m-1}(s-n)|n⟩⟨n|,$$

$${\mathbf{S}}^{+}=\sum _{n=0}^{m-2}\left([k+1{]}_q+\frac{\alpha \beta }{[2s-k{]}_q}\right)|n+1⟩⟨n|+\alpha |m-1⟩⟨0|,$$

$${\mathbf{S}}^{-}=\sum _{n=0}^{m-2}[2s-k{]}_q|n⟩⟨n+1|+\beta |0⟩⟨m-1|,$$

$$Y(\lambda )=[{\partial }_{\beta }T(\lambda ){]}_{\beta =0};T(\lambda )=\mathrm{tr}{}_a(L(\lambda {)}_{a,1}\cdot L_{a,2}(\lambda )\cdot ...\cdot L_{a,n}(\lambda )).$$

$$[\sum _i{s}_i^z,Y(\lambda )]=mY(\lambda ).$$

$$⟨O_k(t)O_l⟩=C_{k,l,m}\cos (\omega t),$$

$$C_{k,l,m}=\left|2\int \mathrm{d}\lambda f_{k,m}(\lambda )⟨Y(\lambda )O_l⟩\right|.$$

$${\kappa }_m(\lambda ,\mu )=⟨Y(\lambda )Y^{†}(\mu )⟩=\frac{\sin \lambda \sin \bar{\mu }}{2{\sin }^2\eta }\frac{\sin (\lambda +\bar{\mu })}{\sin (m(\lambda +\bar{\mu }))}.$$

$$K_m(\lambda ,\mu )={\kappa }_m(\frac{\pi }{2}+\mathrm{i}\lambda ,\frac{\pi }{2}+\mathrm{i}\mu )=\frac{\cosh (\lambda )\cosh (\mu )}{2{\sin }^2\eta }\frac{\sinh (\lambda -\mu )}{\sinh (m(\lambda -\mu ))}.$$

$$O_{k,m}(\lambda )=⟨O_k{Y}^{†}⟩=C(m,k)(\cosh \lambda {)}^{-k+2}.$$

$$\int \mathrm{d}\lambda \frac{\cosh (\lambda )\cosh (\mu )}{2{\sin }^2\eta }\frac{\sinh (\lambda -\mu )}{\sinh (m(\lambda -\mu ))}{f}_{k,m}(\lambda )=C(m,k)(\cosh \mu {)}^{-k+2},$$

$$\int \mathrm{d}\lambda \frac{\sinh (\lambda -\mu )}{\sinh (m(\lambda -\mu ))}{f}_{k,m}^{\prime }(\lambda )=(\cosh \mu {)}^{-k+1},$$

$${\tilde{f}}_{k,m}^{\prime }(\zeta )=\frac{1}{\sqrt{2\pi }}\frac{{\tilde{A}}_k(\zeta )}{{\tilde{K}}_m(\zeta )}.$$

$${\tilde{A}}_{k+1}(\zeta )=\{\begin{array}{ll}{\Pi }_{l=1}^{k/2-1}({\zeta }^2+(2l{)}^2)\zeta \frac{\sqrt{\frac{\pi }{2}}\text{csch}\left(\frac{\pi \zeta }{2}\right)}{(k-1)!};& k\in 2\mathbb{Z}\\ {\Pi }_{l=1}^{(k-1)/2}({\zeta }^2+(2l-1{)}^2)\frac{\sqrt{\frac{\pi }{2}}\text{sech}\left(\frac{\pi \zeta }{2}\right)}{(k-1)!};& k\in 2\mathbb{Z}+1\end{array}$$

$${\tilde{K}}_m(\zeta )=\frac{\sqrt{\frac{\pi }{2}}\tan \left(\frac{\pi }{m}\right)}{m\left(\frac{\cosh \left(\frac{\pi \zeta }{m}\right)}{\cos \left(\frac{\pi }{m}\right)}+1\right)},$$

$$C_{k,l,m}=\left|4C(m,l){\sin }^2\eta \int \frac{\mathrm{d}\lambda }{\cosh (\lambda )}\int \mathrm{d}\zeta \exp (\mathrm{i}\zeta \lambda ){\tilde{f}}_{k,m}^{\prime }(\zeta )\int \mathrm{d}{\zeta }^{\prime }\exp (-\mathrm{i}{\zeta }^{\prime }\lambda ){\tilde{A}}_{l-1}({\zeta }^{\prime })\right|.$$

$$C_{k,l,m}=\left|4\pi C(m,l){\sin }^2\eta \int \mathrm{d}\zeta \mathrm{d}{\zeta }^{\prime }\frac{{\tilde{f}}_{k,m}^{\prime }(\zeta ){\tilde{A}}_{l-1}({\zeta }^{\prime })}{\cosh (\frac{\pi }{2}(\zeta -{\zeta }^{\prime }))}\right|.$$

$$\mathbf{X}(\eta )=\frac{1}{2}\left(\sum _{k=0}^{m-2}\sin (k\eta )|k⟩⟨k+1|+\sin ((k+1)\eta )|k+1⟩⟨k|\right).$$

Full text

Isolated systems consisting of many interacting particles are generally assumed to relax to a stationary equilibrium state whose macroscopic properties are described by the laws of statistical physics. This has been confirmed by a large amount of recent theoretical (e.g. Eisert et al. (2015); Sotiriadis and Martelloni (2016); Collura et al. (2013); Sotiriadis and Calabrese (2014); Calabrese et al. (2016, 2011); Manmana et al. (2007); Rossini et al. (2009); Caux and Konik (2012); Rigol (2014); Mestyán et al. (2017); Polkovnikov et al. (2011); Gogolin and Eisert (2016); Eckstein et al. (2009); Alba and Calabrese (2017); Calabrese et al. (2012); Rigol et al. (2008)) and experimental (e.g. Cheneau et al. (2012); Langen et al. (2013); Gring et al. (2012); Chen et al. (2011); Meinert et al. (2013); Langen et al. (2015); Clos et al. (2016)) work, particularly focusing on the Heisenberg spin chain in quench setups (e.g. Wouters et al. (2014); Pozsgay et al. (2014); Bertini et al. (2016); Mestyán et al. (2015); Piroli et al. (2016); Ilievski et al. (2015a); Collura et al. (2015); Ilievski et al. (2017)). On the other hand time crystals describe a phase of mater which never relaxes to stationarity and breaks the continuous time-translation symmetry (TTS) in analogy with the continuous space translation symmetry breaking in ordinary crystals. Historically the research into quantum time crystals was instigated by an intriguing possibility that a system at zero temperature could exhibit perpetual motion Wilczek (2012), however this has subsequently been disputed Bruno (2013a, b); Watanabe and Oshikawa (2015), leading to possible generalisations to finite temperature, or far from thermal equilibrium. Despite the large amount of work on Floquet (breaking the discrete time translation symmetry) or dissipation induced time crystals Else et al. (2016); Lazarides et al. (2014); Choi et al. (2017); Bordia et al. (2017); Sacha and Zakrzewski (2017); Yu et al. (2019); Zhu et al. (2019); Gambetta et al. (2019); Giergiel et al. (2019); Gambetta et al. (2019); Cosme et al. (2019); Buča et al. (2019); Iemini et al. (2018); Dogra et al. (2019); Barberena et al. (2019); Lledó et al. (2019) such behavior was believed to be impossible to realize in isolated many-body systems Bruno (2013a, b); Watanabe and Oshikawa (2015); Volovik (2013), which typically relax to stationary states depending only on few parameters, such as energy and particle number D’Alessio et al. (2016); Vidmar and Rigol (2016).

Despite numerous studies on relaxation in many-body quantum systems, there have been no results on spontaneous time translation symmetry breaking to date, except for mean field models (e.g. in the non-interacting limit Syrwid et al. (2017), fully connected models Sciolla and Biroli (2010)), or the case of spin precession. These systems are not believed to have the required eigenfrequency complexity for local observables to relax to stationarity D’Alessio et al. (2016); Essler and Fagotti (2016). By contrast, we will focus on strongly correlated systems where one does a priori expect relaxation to stationarity Eisert et al. (2015); Sotiriadis and Martelloni (2016); Collura et al. (2013); Sotiriadis and Calabrese (2014); Calabrese et al. (2016, 2011); Manmana et al. (2007); Rossini et al. (2009); Caux and Konik (2012); Rigol (2014); Mestyán et al. (2017); Polkovnikov et al. (2011); Gogolin and Eisert (2016); Eckstein et al. (2009); Alba and Calabrese (2017); Calabrese et al. (2012); Rigol et al. (2008). The absence of such results might be expected in the light of the no-go theorem Watanabe and Oshikawa (2015), however there is a crucial defining property underlying its derivation, which we relax. It assumes that the system should exhibit long range spatial correlations, which are not relevant to time-translation symmetry breaking, or non-stationarity.

In this letter we show that systems can indeed fail to relax, and relate this type of behavior to extensive dynamical symmetries, which are local in space and have a periodic dependence on time. As a consequence of dynamical symmetries systems support time-dependent statistical ensembles, and fail to equilibrate after a quantum quench. Non-stationarity shows up also on the level of dynamical response functions which describe the behavior near equilibrium, and its stability under small perturbations.

We demonstrate the effects of local dynamical symmetries for the one-dimensional Heisenberg spin chain and study the stability under the integrability breaking perturbation. This paradigmatic model is used to describe many experimentally relevant situations including organic componounds Blundell and Pratt (2004), various materials Breunig et al. (2013), cold atom implementations Gross and Bloch (2017), and quantum dots Zwanenburg et al. (2013).

Quantum time crystals. Watanabe and Oshikawa defined quantum time crystals as an interacting system exhibiting persistent oscillations, which can be probed by the auto-correlation function $\tilde{f}(t)=\frac{1}{V^{2}}\langle O(t)O\rangle$, where $V$ is the volume of the system and $O$ an extensive observable Watanabe and Oshikawa (2015). They consider this auto-correlation function as a perturbation to equilibrium and show that it is time independent at zero and finite temperature. Their definition might not capture physically measurable persistent oscillations because the time dependent (connected) part of the function $\tilde{f}(t)$ vanishes in the absence of long range correlations in large volume limit at all times. We instead propose a definition where the connected auto-correlation function is initially normalized $f(t)=\frac{1}{\langle O^{2}\rangle}\langle O(t)O\rangle$, which probes TTS breaking if equilibrium ensemble is perturbed, even in the absence of long range correlations. Following such perturbation quantum time crystals will never reach stationarity.

While there are different ways of identifying the many-body nature of the phenomenon, the definition we use here is that single-body observables relax to stationarity, while some of the many-body observables oscillate persistently. This way the oscillation of the higher point correlation functions cannot be attributed to the single-body oscillations.

Recently, most of the work on quantum time-crystals has focused on the discrete time translation symmetry breaking in Floquet systems Else et al. (2016); Lazarides et al. (2014); Choi et al. (2017); Bordia et al. (2017); Sacha and Zakrzewski (2017); Yu et al. (2019); Zhu et al. (2019); Gambetta et al. (2019); Giergiel et al. (2019); Gambetta et al. (2019). In this case TTS breaking is studied in the out-of-equilibrium quench setup. In this setup the system is prepared in a generic pure state $|\psi\rangle$ and then allowed to evolve under the action of a Hamiltonian. The system is identified as a discrete quantum time crystal if the dynamics breaks the discrete time symmetry of the driving period $T$ with a subharmonic response $\langle\psi|o(t+nT)|\psi\rangle=\langle\psi|o(t)|\psi\rangle$, for some integer $n>1$, and $\langle\psi|o(t+t_{1})|\psi\rangle\neq\langle\psi|o(t)|\psi\rangle$ for $t_{1}<nT$ Else et al. (2016). We make the analogous identification for continuous time evolution by requiring that $\langle\psi|o(t+T)|\psi\rangle=\langle\psi|o(t)|\psi\rangle$ and $\langle\psi|o(t+t_{1})|\psi\rangle\neq\langle\psi|o(t)|\psi\rangle$, for $t_{1}<T$ and for some $T$.

Extensive dynamical symmetries. An important insight into phenomena of equilibration is provided by the eigenstate thermalization hypothesis in generic systems Deutsch (1991); Srednicki (1994); Rigol et al. (2008); D’Alessio et al. (2016), or generalized eigenstate thermalization hypothesis (GETH) in integrable systems Cassidy et al. (2011). It states that offdiagonal elements of local observables in eigenbasis of local Hamiltonian vanish exponentially in thermodynamic limit, and that their expectation values in a given eigenstate are smooth functions of conserved quantities. Assuming the validity of GETH the system is expected to locally relax to the maximal entropy, or generalized Gibbs ensemble (GGE) $\rho_{GGE}=\exp(-\sum_{j}\mu_{j}X_{j})$, following the quantum quench. The set of chemical potentials $\mu_{j}$ is obtained by matching the expectation values of extensive conservation laws $X_{j}$ Ilievski et al. (2016a) in the ensemble and the initial state Essler and Fagotti (2016). The situation is very different if the system possesses an additional set of extensive dynamical symmetries $Y$ and $Y^{\dagger}$. Such quantities satisfy a simple closure (or eigenoperator) condition

[TABLE]

which leads to periodic evolution $Y(t)=\exp({\rm i}\omega t)Y(0)$ and $Y(t)^{\dagger}=\exp(-{\rm i}\omega t)Y(0)^{\dagger}$. In any isolated system there is a large number of operators satisfying condition (1), however in general they are highly nonlocal, and as such have no effect on local physics on large timescales. For simplicity we will discuss the case with a single frequency $\omega$, but the generalization to multiple frequencies is straightforward.

Heisenberg model and extensive dynamical symmetries. Throughout the letter we will consider the anisotropic Heisenberg Hamiltonian

[TABLE]

as an example of the quantum time crystal. In equation (LABEL:Heisenberg) we introduced the spin$-{\textstyle\frac{1}{2}}$ operators $s^{x,y,z}$, anisotropy $\Delta$, hopping amplitude $J$, the magnetic field $h$, and the integrability breaking parameter $\alpha$, which is set to [math] except when otherwise specified. One of the crucial aspects of the Heisenberg model, which has a paramount effect on physical properties are the extensive conservation laws Prosen ; Zadnik et al. (2016); Ilievski et al. (2015b); Prosen and Ilievski (2013). These conservation laws are extensive operators that commute with the Hamiltonian and have an overlap with physical local quantities. Their effects range from the absence of thermalization to the ideal energy and spin conductivity at any temperature. Despite the absence of thermalization, the Heisenberg model has in recent years served as a testbed for studying equilibration properties of strongly interacting systems Wouters et al. (2014); Pozsgay et al. (2014); Bertini et al. (2016); Mestyán et al. (2015); Piroli et al. (2016); Ilievski et al. (2015a); Collura et al. (2015); Ilievski et al. (2017). In what follows we will show that in the easy plane regime $-1<\Delta<1$, it, in general, never reaches equilibrium if $h\neq 0$.

This can be seen as a consequence of semi-cyclic quantities, which were introduced in Zadnik et al. (2016) (see also SM ). While they commute with Hamiltonian (LABEL:Heisenberg) in the absence of the field $h=0$, we will show SM that they satisfy the closure condition (1) for any value of the field $h$. Interestingly, the frequency of their oscillations $\omega=hm$ is a discontinuous function of the anisotropy parameter $\Delta=\cos(\frac{\pi n}{m})$, with $n\in 2\mathbb{N},\ m\in 2\mathbb{N}+1$ (see FIG. 1). Moreover the explicit structure of dynamical symmetries depends finely on the exact value of anisotropy $\Delta$, since they are comprised of densities which have a surplus of $m$ local operators $s^{+}$. Conversely, quantities $Y^{\dagger}$ have a surplus of local operators $s^{-}$. This means that, for instance, at $\Delta=-\frac{1}{2}$ we will observe the persistent oscillations of the three point transverse correlation function $s^{x}_{1}s^{x}_{2}s^{x}_{3}$, while the oscillations of this observable in the infinite time limit will be absent at any other value of $\Delta$. Physically, such observables correspond to correlations of the three-site measurement statistics - the average measured value of each individual spin relaxes according to standard statistical physics, but the measured values will be such that on average their product oscillates in time. Alternatively, they may be thought of as oscillations of the higher moments of the $m$-site quantum fluctuations (e.g. $\langle(s^{x}_{1}+s^{x}_{2}+s^{x}_{3})^{3}(t)\rangle$). Note that the fact that quantities responsible for oscillations Y do not exist at the non-interacting point of the model with $\Delta=0$ solidifies the argument that the oscillations are a genuine many-body phenomena.

Time dependent generalized Gibbs ensemble (tGGE). In order to understand how the dynamical symmetries affect late time dynamics of local observables after the quantum quench, we consider discrete time dynamics induced by a Hamiltonian $H$ with a period $\frac{2\pi}{\omega}$,

[TABLE]

which renders $Y$ and $Y^{\dagger}$ conserved. The stationary maximum entropy ensembles for the stroboscopic dynamics can be obtained from the entropy maximization procedure, which has to respect the conservation of all conserved quantities $X_{j}$, as well as dynamical symmetries $Y$ and $Y^{\dagger}$. This leads to the GGE description $\rho_{GGE}=\exp(-\sum_{j}\mu_{j}X_{j}-\mu_{Y}Y-\bar{\mu}_{Y}Y^{\dagger})$. If conserved quantities do not commute this might in principle require redefinition of ensembles Doyon (2017); Fagotti (2014). Stroboscopic time evolution of the state $|\psi(t)\rangle$ leads to different maximum entropy states for $t\in[0,\frac{2\pi}{\omega})$, which take the form of time dependent generalized Gibbs ensemble (tGGE)

[TABLE]

The values of the chemical potentials $\mu_{j}$, $\mu_{Y}(t)$ can be fixed in the following way. The maximum entropy non-stationary ensemble which correctly reproduces the initial value of $X_{k}$, and the dynamics of dynamical symmetries $Y_{k}$ is obtained by requiring

[TABLE]

An important observation is that if the following set of relations holds

[TABLE]

the constraints (5) can be satisfied by the time dependent chemical potentials $\mu_{Y}(t)=\mu_{Y}(0)\exp(-{\rm i}\omega t)$.

Dynamical response functions in thermal equilibrium. Here we focus on the response of the system in equilibrium. Extensive dynamical symmetries have profound consequences for the asymptotic behavior of dynamical susceptibilities as well. In the large time limit we expect that a local observable $O(t)$ can be represented as a linear combination of conserved quantities and dynamical symmetries

[TABLE]

for the sake of calculation of the dynamical susceptibilities $\langle O_{1}(t)O_{2}\rangle$. This equality is valid only on the hydrodynamic level and in the long-time limit. Here we restricted the discussion to the case of $\langle O_{1}\rangle=\langle O_{2}\rangle=0$, which can be relaxed by considering connected correlation functions. In general the coefficients $\alpha_{Y}^{O}$ and $\alpha_{j}^{O}$ depend, not only on observable $O$, but also on the thermal ensemble $\langle\bullet\rangle$. If these coefficients vanish, the observable is not expected to oscillate.

In order to obtain the dynamics of temporal correlation functions in Heisenberg model $\langle O_{k}(t)O_{l}\rangle$ with $O_{k}=\sum_{i}s^{x}_{i}...s^{x}_{i+k-1}$, we will use the ansatz (7), specializing to the infinite temperature ensemble $\langle\bullet\rangle=\frac{\,{\rm tr}\,(\bullet)}{\,{\rm tr}\,(\mathbbm{1})}$. Conversely, the results provide an asymptotic solution of the quench protocol for any initial state of the form $\rho=\sum_{k\in 2\mathbb{Z}+1}a_{k}\,O_{k}$. A particular choice of $O_{k}$ was made due to the non-zero overlap with the Y quantities of the Heisenberg XXZ spin chain Zadnik et al. (2016) in the sense of (7).

A starting point of the derivation is expression (7), where we take into account that the set of extensive conserved quantities and dynamical symmetries in the Heisenberg model is comprised of three single parameter families: unitary conservation laws $X_{s}(\lambda)$, non-unitary conservation laws $Z(\lambda)$, and semi-cyclic dynamical symmetries $Y(\lambda)$ with the frequency which is independent of $\lambda$ Ilievski et al. (2016a). In canonical ensemble the observable $O_{k}$ can be described asymptotically as

[TABLE]

Note that the judicious choice of the $O_{k}$ and the canonical ensemble means that the charges $X_{k}(\lambda)$ and $Z(\lambda)$ are not relevant via (7). The function $f(\lambda)$ and its complex conjugate $\bar{f}(\lambda)$ can be obtained by projecting the expression (8) onto the set of dynamical symmetries

[TABLE]

Specializing to the infinite temperature, the overlaps $\langle O_{k}Y^{\dagger}(\lambda^{\prime})\rangle$, and kernels $\langle Y(\lambda)Y^{\dagger}(\lambda^{\prime})\rangle$ were obtained in Zadnik et al. (2016). The system of equations (9) can be reduced to the convolution equation and solved in Fourier space (see SM for more details). We also note that this choice of temperature implies that $\langle O^{2}\rangle\propto V$. Using the solution $f(\lambda)$, we can calculate the susceptibility matrix $\langle O_{k}(t)O_{l}\rangle\underset{t\to\infty}{=}C\cos(\omega t)$, with $C=\left|2\int{\rm d}\lambda f(\lambda)\langle Y(\lambda)O_{l}\rangle\right|$. For instance, setting the anisotropy parameter to $\Delta=-\frac{1}{2}$, and $k=l=3$ we get $C=\frac{1}{64}\left(\frac{27\sqrt{3}}{\pi}-8\right)$. In general the result is expressed as a double integral of elementary functions (see SM ), and the constant which can be calculated efficiently.

Numerical results. In Fig. 2 a) we plot the time evolution of the autocorrelation function $C(t)=\frac{\langle O(t)O\rangle}{\langle O^{2}\rangle}$, for the observable $O=\sum_{j}s_{j}^{x}s_{j+1}^{x}s_{j+2}^{x}$, and the infinite temperature ensemble $\langle\bullet\rangle=\frac{\,{\rm tr}\,(\bullet)}{\,{\rm tr}\,(\mathbbm{1})}$. Clearly $C(t)$ does not equilibrate, and we can compare the numerical result with the analytical calculation $C(t)\underset{t\to\infty}{=}\frac{1}{64}\left(\frac{27\sqrt{3}}{\pi}-8\right)\cos(3ht)$, obtaining a perfect agreement. As predicted by theory, we observe no oscilations of the transverse magnetization $O=\sum_{j}s_{j}^{x}$.

In Fig. 2 b) we plot the time dependence of a local three point correlation function $o=s_{j}^{+}s_{j+1}^{+}s_{j+2}^{+}+s_{j}^{-}s_{j+1}^{-}s_{j+2}^{-}$ from the ferromagnetic initial state maximally polarized in the $x$ direction,

[TABLE]

as well as the one-point function $o=s_{j}^{+}+s_{j}^{-}$ (single body observable). It shows the relaxation of the one-point function and persistent oscillations in the three point function, illustrating the many body nature of the time crystalline behaviour.

As shown in Fig. 2 c) we find numerically that the oscillations of relevant observables persist for a long time with an altered amplitude, following small to intermediate perturbations of the anisotropy. At $\Delta\neq-0.5$ the three-point operator $s^{x}_{j}s^{x}_{j+1}s^{x}_{j+2}$ no longer has an overlap with the $Y$ ($\,{\rm tr}\,(Ys^{x}_{j}s^{x}_{j+1}s^{x}_{j+2})=0$, ) and thus does not pertain to any dynamical symmetry (see Eqs. (5),(7)). Only at a relatively large value of perturbation $\delta\Delta=0.2$ do we see significant damping of the amplitude on the accessible time scales. The oscillations are stable also with respect to the integrability breaking term $\alpha$ on accessible timescales as shown in Fig. 2 d). We use the general results of Kollar et al. (2011); Mori et al. (2018) to find that the $Y$ operators are conserved under stroboscopic time evolution in Eq. (3) up to at-least second-order in perturbation strength (see also Mallayya et al. (2019)). This implies that the oscillations decay no faster than $\exp(-tJ{\cal O}(\alpha^{2}))$. Though this may possibly be improved as in the cases where the bounds on prethermalization are exponential Abanin et al. (2017); Mori et al. (2018); Brandino et al. (2015).

Experimental realization. Due to the demonstrated stability, we expect that in current quantum cold-atom simulations of the XXZ spin chain, such as those done in I. Bloch’s group Fukuhara et al. (2013), the lattice depth can be sufficiently tuned to make the dynamics fast enough compared to integrablity-breaking effects to observe oscillations.

Measurement of local on-site equal-time many-body correlation functions, such as the ones we study, is available through quantum gas microscopes for cold atom systems Gross and Bloch (2017). For experiments an important discovery is that oscillations can be observed for a quench from the ferromagnetic initial state (10), which can be engineered Gross and Bloch (2017). Preparation and measurement of auto-correlation functions is more involved, but can be achieved through the use of ancilla qubits in Rydberg atoms (see SM for details). Our results could potentially also have far-reaching applications in quantum metrology Giovannetti et al. (2011), as they offer, in theory, infinitely sensitive and precise measurement of the system anisotropy. In cold atom simulations this can be directly related to the strength of the external magnetic field used to achieve Feshbach resonance of the spin-spin interaction Gross and Bloch (2017). In this regard, an important observation is that the amplitude seems to be less effected by integrability breaking, then by the change of anisotropy.

Conclusion: Numerous questions remain open. While we have addressed the question of stability to perturbations from a practical perspective, stability to all orders remains an open problem, related to the long standing question of the existence of KAM theorem in systems with infinitely many degrees of freedom Kappeler and Pöschel (2013). That being said, the crucial ingredient for oscillations is not integrability itself, but rather local or extensive quantities satisfying the relation (1). Importantly, we were able to identify similar quantities in topological models Wen (2003), however we postpone the in-depth discussion on this topic to later publications. The glimmers of similar dynamical symmetries have also been identified in locally constrained models exhibiting quantum many-body scars, preventing the systems from relaxing for certain initial conditions Turner et al. (2018); Choi et al. (2019); Bernien et al. (2017).

Some questions remain also from the standpoint of integrable systems. Here we only focused on the lowest frequency of oscillations at a given $\Delta$, while in general the state $\rho_{tGGE}$ should support a complete harmonic spectrum $\omega=khm$, for $k\in\mathbb{N}$. Furthermore, due to the noncommutativity of the conserved quantities and dynamical symmetries, subtleties might arise in obtaining the correct form of $\rho_{tGGE}$ Doyon (2017). The answers to these questions should be attainable by extending thermodynamic Bethe ansatz description Ilievski et al. (2016b), to include additional quantities $Y$. The existence of these quantities implies that the standard GGE description is in fact incomplete even in the absence of persistent oscillations ($h=0$). Another exciting question is whether the dynamical symmetries, and absence of many-body equilibration has a counterpart in the realm of classical physics. Otherwise, the phenomenon would constitute one of the first many-body quantum effects that can be observed in macroscopic systems on large space and time scales solely due to the extensive many-body nature of the $Y$ operators.

S1 Semi-cyclic dynamical symmetries

A starting point for constructing the conservation laws in Heisenberg model is the Lax matrix ChargesReview_ap

[TABLE]

where the operators acting on the auxiliary space ${\bf S}^{\alpha}$ have to satisfy the $U_{q}(\mathfrak{sl}_{2})$ quantum group relations QGBook_ap

[TABLE]

The constant $\eta$ relates to the anisotropy as $\Delta=\cos(\eta)$. There are three distinct families of extensive conserved quantities in the anizotropic Heisenberg model if the magnetic field is absent $h=0$. For our purpose, the relevant family is the family of semicyclical conservation laws YCharges_ap , which exist for the value of $\eta=\frac{\pi n}{m}$, and cover the interval $-1<\Delta<1$ densely. These conservation laws are obtained from the (semi)cyclic representation of $U_{q}(\mathfrak{sl}_{2})$ korf

[TABLE]

where the $q$-deformation of a number $[a]_{q}$ corresponds to $[a]_{q}=\frac{q^{a}-q^{-a}}{q-q^{-1}}$. Note that quantities are conserved only for even $n$ and odd $m$. Conserved quantities are extensive in the region ${\,{\rm Re}\,}(\lambda)\in\{\frac{\pi}{2}-\frac{\pi}{m},\frac{\pi}{2}+\frac{\pi}{m}\}$, if the representation parameters $s=0$, and the derivative of the transfer matrix $T(\lambda)$ is taken with respect to parameter $\beta$ at $\beta=\alpha=0$ YCharges_ap

[TABLE]

What is important for our purpose is that this set of conservation laws does not commute with magnetization and as a consequence with Heisenberg Hamiltonian if $h\neq 0$. It is, however, rather simple to show that they satisfy the closure condition from the main text (1). This is a consequence of operators $Y(\lambda)$ having a surplus of operators $s^{+}$. The local densities of magnetization are supported on a single site, implying that we can calculate the commutator by considering a single site commutators of $s^{z}$ with $s^{\pm,z,0}$, $[s^{z},s^{\pm,z,0}]=\alpha^{\pm,z,0}s^{\pm,z,0}$. Since $\alpha^{z,0}=0$, and $\alpha^{\pm}=(-1)^{\frac{1\mp 1}{2}}$ and $Y$ comprises operators with a fixed surplus $m$ of $s^{+}$ we have

[TABLE]

Let us now briefly discuss some additional properties of quantities $Y(\lambda)$. First of all, it is obvious that $Y(\lambda)Y^{\dagger}(\lambda)$ is conserved. Interestingly enough, for finite sizes quantities $Y(\lambda)$ do not commute with $Y^{\dagger}(\lambda)$. This property is interesting, since the commutator of two local operator is a local operator, and we might expect that the commutator of two quasilocal operators is quasilocal as well. This, however, is nontrivial to show. Furthermore, the commutator of the operators $A(\lambda,\mu)=[Y(\lambda),Y^{\dagger}(\mu)]$, is spin flip $S=\sigma^{x}_{1}\otimes\sigma^{x}_{2}\otimes...\otimes\sigma^{x}_{n}$ antisymmetric $SA(\lambda,\mu)S=-A(\lambda,\mu)$ by construction, and conserved for any value of the field $h$. Provided that the commutator is extensive it should be expressible in terms of charges $Z(\lambda)$ Prosen_periodic , $A(\lambda,\mu)=\int{\rm d}\nu\,h(\nu)Z(\nu)$, since the set of the charges $Z(\lambda)$ seems to be complete Zchargi_app ; Prosen_periodic ; EnejJacopo ; Delucacharges ; LZP . Furthermore, one can check that quantities $Y(\lambda)$ do not commute with $A(\lambda,\mu)$ for small systems, and could in principle lead to new $U(1)$ symmetry breaking charges. We, however, do not think that this is the case for two reasons: firstly numerically obtained dynamical susceptibility describes the dynamical correlation functions perfectly for $\Delta=-\frac{1}{2}$, and secondly semicyclic charges are closely related to the current carrying charges $Z(\lambda)$, for which, as already mentioned, there is a lot of evidence indicating their completeness.

S2 Time dependent dynamical susceptibilities

In this section we discuss how to analyticaly obtain the asymptotic values of time dependent susceptibilities (autocorrelation functions) for observables $O_{k}=\sum_{i}\sigma^{x}_{i}...\sigma^{x}_{i+k-1}$, using the method of hydrodynamical projection hydro . First we focus on the infinite temperature case, and later discuss the general GGE case. The starting point is the expression (8) from the main text, which after the multiplication with $O_{l}$ from the right yields

[TABLE]

up to the phase factor, and $m$ corresponds to the parametrization of anisotropy $\Delta=\cos(\tfrac{\pi n}{m})$. The constant $C_{k,l,m}$ is

[TABLE]

Note that we reintroduced indices $k$ and $m$ in the function $f$, and that $Y(\lambda)$ are $m$ dependent as well.

Let’s now focus on solving the equation (9) from the main text, and in turn provide an explicit expression for constants $C_{k,l,m}$. The starting point is the kernel derived in YCharges_ap

[TABLE]

Using the conjecture that the complete set of charges can be obtained as a power expansion of the continuous family $Y(\lambda)$ along the line $\lambda=\frac{\pi}{2}+{\rm i}\,t$ wrt. $t$, we can restrict the discussion to the family of quantities $Y(\lambda)$ with a single real parameter

[TABLE]

The overlap of $Y(\lambda)$ with observable $O_{k}$ reads YCharges_ap

[TABLE]

This leads to a set of Fredholm equation

[TABLE]

which can be reduced to the set of convolution equations by rescaling the function $f_{k,m}(\lambda)=C(m,k)\frac{2\sin^{2}\eta}{\cosh(\lambda)}f_{k,m}^{\prime}(\lambda)$:

[TABLE]

and solved in Fourier space

[TABLE]

Here we introduced the Fourier transform of the rhs (S15) $\tilde{A}_{k}(\zeta)=FT[(\cosh\lambda)^{-k+1}]$, and conjecture its analytical form

[TABLE]

The kernel in Fourier space reads $\tilde{K}_{m}(\zeta)=FT[\frac{\sinh(\lambda)}{\sinh(m\lambda)}]$,

[TABLE]

This implies that we can write the dynamical structure constants as

[TABLE]

and finally the integration over ${\rm d}\lambda$ yields

[TABLE]

Specializing to the value $\eta=\frac{2\pi}{3}$, we conjecture the value of the constant in the expression for $O_{k,m}(\lambda)$ is $C(3,k)=\sqrt{3}\,3^{k-3}/16^{2^{(k-m)}}$. In particular, if we set $k=l=3$, we get the final result $C_{3,3,3}=\frac{1}{64}\left(\frac{27\sqrt{3}}{\pi}-8\right)$. The function $C(k,m)$ can be obtained by evaluating $k$-th power $C(k,m)=(-1)^{(k-m)}\langle m-1|\mathbf{X}(\eta)^{k}|1\rangle$ of matrix

[TABLE]

In general we can show for $O$ with the smallest support with non-zero constant (namely, $k=m$) that $C(m,m)=\frac{1}{2^{m}}\prod_{\mu=2}^{m-1}\sin(\mu\eta)$.

To illustrate another value of $\Delta$ in Fig. S3 we show the autocorrelation function for $O=\sum_{j}s^{x}_{j}s^{x}_{j+1}s^{x}_{j+2}s^{x}_{j+3}s^{x}_{j+4}$ at $\Delta=\cos\left(\tfrac{2\pi}{5}\right)$.

S3 Measurement of correlation functions

In this discussion we focus on cold atom experimental setups coldexp_ap , though other implementations may be possible. The measurement of the local time correlation functions, such as $o=\langle s^{+}_{i}s^{+}_{i+1}s^{+}_{i+2}\rangle+h.c$ from the main text, can be achieved by on-site quantum gas microscopy, as is discussed for two-point correlation functions in coldexp_cor in the Fermi-Hubbard model. Combining existing techniques for simulating the XXZ spin chain with cold atoms coldexp_ap with quantum gas microscopy would be required. Measurement of the examples for $o$ from the main text can be done by collecting the statistics of repeated measurements from the same initial state by measuring the local Hermitian observables making up $o$, e.g. $s^{x}_{i}s^{y}_{i}s^{x}_{i}$, etc.

Measurement of the autocorrelation functions of the form $\langle O(t)O\rangle/{\langle O\rangle^{2}}$ is more involved, but it could be achieved in three ways utilizing current technologies.

Firstly, one may collect the statistics of the joint probability distribution for $O$ at different times. Preparing the initial infinite temperature state could be accurately achieved by preparing a very high initial temperature state. This should be followed by an application of the external magnetic field in the z-direction to avoid introducing $s^{z}_{i}$ terms in the initial density matrix.

Secondly, we may also easily observe that measuring $\langle O(t)O\rangle$ is equivalent to measuring $\langle O(t)\rangle$ for a quench from an initial density matrix of the form $\rho(0)=\exp{\left(\mu O\right)}$ for very small $\mu$ (linear response regime). Such an initial state could be prepared with Rydberg atoms. Introducing ancilla qubits to the XXZ spin chain setup would mimic local three-site interactions of the form $O=\sum_{i}s^{x}_{i}s^{x}_{i+1}s^{x}_{i+2}$, as was proposed for a transverse-field Ising model, but in 2D Alex . By tuning the external field as discussed in Alex we can make the $a/\mu$ in $H=aH_{\rm XXZ}+\mu O$ very small for small $\mu$, thus realizing a state close to $\rho(0)$. The field would then be quenched Alex such that $a$ becomes large and a homogenous magnetic field in the z-direction introduced. One may then measure $O(t)$ with on-site techniques, as discussed previously.

Thirdly, one can use the method of single-qubit interferometry inter to measure the autocorrelation by measuring a probe qubit as discussed in Juha for measuring the Green’s function in a strongly-correlated electronic system. This would require adapting the same protocol for spins.