Emergent prethermalization signatures in out-of-time ordered correlations

TL;DR

This paper reports on nuclear magnetic resonance experiments that observe prethermalization signatures in out-of-time ordered correlations within a spin chain, revealing emergent conserved quantities and advancing understanding of quantum thermalization.

Contribution

The study introduces a novel protocol for measuring out-of-time ordered correlations and demonstrates the emergence of prethermalization in a controllable quantum system.

Findings

Emergence of conserved quantities at large magnetic fields.

Early saturation of correlations indicating prethermalization.

Validation of Hamiltonian engineering techniques for probing quantum dynamics.

Abstract

How a many-body quantum system thermalizes --or fails to do so-- under its own interaction is a fundamental yet elusive concept. Here we demonstrate nuclear magnetic resonance observation of the emergence of prethermalization by measuring out-of-time ordered correlations. We exploit Hamiltonian engineering techniques to tune the strength of spin-spin interactions and of a transverse magnetic field in a spin chain system, as well as to invert the Hamiltonian sign to reveal out-of-time ordered correlations. At large fields, we observe an emergent conserved quantity due to prethermalization, which can be revealed by an early saturation of correlations. Our experiment not only demonstrates a new protocol to measure out-of-time ordered correlations, but also provides new insights in the study of quantum thermodynamics.

Figures

Equations

$$H_{\text{TDip}}=u\sum _{j,k>j}{J}_{jk}\left[S_y^j{S}_y^k-\frac{1}{2}(S_x^j{S}_x^k+S_z^j{S}_z^k)\right]+g\sum _j{S}_z^j$$

$$S_{\varphi }=2^{-L}\text{Tr}[e^{-i\varphi {\mathcal{O}}_{\mathbf{n}}}\delta \rho (t)e^{i\varphi {\mathcal{O}}_{\mathbf{n}}}\delta \rho (t)].$$

$$C_{{\mathcal{O}}_{{\mathbf{n}}^{\prime }}{\mathcal{O}}_{\mathbf{n}}}(t)=\frac{4}{L}⟨|[{\mathcal{O}}_{{\mathbf{n}}^{\prime }}(t),{\mathcal{O}}_{\mathbf{n}}]{|}^2{⟩}_{\beta =0}=\sum _q{q}^2{I}_q(t)$$

$$H_{\mathrm{pre}}=Z^{\prime }+\sum _{n=1}^{n_M}{ϵ}^n{h}^{(n)},$$

$$Z(Jt=10^3)=\sqrt{2^{L-2}L}\sum _{k=1}^L\sum _{s=1}^{{\zeta }_k}{b}_k^s(Jt){\mathcal{B}}_k^s,$$

$$H_{tot}={\omega }_F\sum _k{S}_z^k+{\omega }_P\sum _{\kappa }{s}_z^{\kappa }+H_F+H_P+H_{FP}$$

$$H_{dip}=\sum _{j<k}\frac{\hslash {\gamma }_j{\gamma }_k}{|{\vec{r}}_{jk}{|}^3}\left[{\vec{S}}_j\cdot {\vec{S}}_k-\frac{3{\vec{S}}_j\cdot {\vec{r}}_{jk}{\vec{S}}_k\cdot {\vec{r}}_{jk}}{|{\vec{r}}_{jk}{|}^2}\right],$$

$$H_F=\frac{1}{2}\sum _{j<k}{J}_{jk}^F(2S_z^j{S}_z^k-S_x^j{S}_x^k-S_y^j{S}_y^k)H_P=\frac{1}{2}\sum _{\lambda <\kappa }{J}_{\kappa \lambda }^P(2s_z^{\lambda }{s}_z^{\kappa }-s_x^{\lambda }{s}_x^{\kappa }-s_y^{\lambda }{s}_y^{\kappa })$$

$$H_{FP}=\sum _{k,\kappa }{J}_{k,\kappa }^{FP}{S}_z^k{s}_z^{\kappa },$$

$$H^{(0)}=\frac{1}{t_c}{\int }_0^{t_c}H(t)dt,H^{(1)}=\frac{-i}{2t_c}{\int }_0^{t_c}d{t}_2{\int }_0^{t_2}d{t}_1[H(t_2),H(t_1)].$$

$$P({\tau }_1,\mathbf{x},{\tau }_2,\mathbf{y},2{\tau }_1,\mathbf{y},{\tau }_2,\mathbf{x},{\tau }_1)P({\tau }_1,\mathbf{x},{\tau }_2,\mathbf{y},2{\tau }_1,\mathbf{y},{\tau }_2,\mathbf{x},{\tau }_1)P({\tau }_1,\bar{\mathbf{x}},{\tau }_2,\bar{\mathbf{y}},2{\tau }_1,\bar{\mathbf{y}},{\tau }_2,\bar{\mathbf{x}},{\tau }_1)P({\tau }_1,\bar{\mathbf{x}},{\tau }_2,\bar{\mathbf{y}},2{\tau }_1,\bar{\mathbf{y}},{\tau }_2,\bar{\mathbf{x}},{\tau }_1)$$

$$P({\tau }_3,\mathbf{y},{\tau }_3,\mathbf{x},2{\tau }_4,\mathbf{x},{\tau }_3,\mathbf{y},{\tau }_3)P({\tau }_3,\mathbf{y},{\tau }_3,\mathbf{x},2{\tau }_4,\mathbf{x},{\tau }_3,\mathbf{y},{\tau }_3)P({\tau }_3,\bar{\mathbf{y}},{\tau }_3,\bar{\mathbf{x}},2{\tau }_4,\bar{\mathbf{x}},{\tau }_3,\bar{\mathbf{y}},{\tau }_3)P({\tau }_3,\bar{\mathbf{y}},{\tau }_3,\bar{\mathbf{x}},2{\tau }_4,\bar{\mathbf{x}},{\tau }_3,\bar{\mathbf{y}},{\tau }_3)$$

$${\tau }_1=\tau (1-u),{\tau }_2=\tau (1+2u),{\tau }_3=\tau (1+u),{\tau }_4=\tau (1-2u)$$

$$H^{(0)}=\pm u\sum _{j<k}{J}_{jk}\left[S_y^j{S}_y^k-\frac{1}{2}(S_x^j{S}_x^k+S_z^j{S}_z^k)\right],$$

$$U_1=e^{-i{H}^{(0)}{t}_c},U_2=e^{-i\varphi Z}{e}^{-i{H}_{\mathrm{Flq}}{t}_c}{e}^{i\varphi Z},U_3=e^{-2i\varphi Z}{e}^{-i{H}_{\mathrm{Flq}}{t}_c}{e}^{2i\varphi Z},\cdots U_n=e^{-i(n-1)\varphi Z}{e}^{-i{H}_{\mathrm{Flq}}{t}_c}{e}^{i(n-1)\varphi Z}$$

$$U(n{t}_c)$$

$$\mathcal{O}=\sum _{q=-L}^{q=L}{\mathcal{O}}_q$$

$$e^{-i\theta P}{\mathcal{O}}_q{e}^{i\theta P}=e^{-iq\theta }{\mathcal{O}}_q,[P,{\mathcal{O}}_q]=q{\mathcal{O}}_q.$$

$${\mathcal{O}}_q=\frac{1}{2\pi }{\int }_0^{2\pi }d\theta e^{iq\theta }{e}^{-i\theta P}\mathcal{O}{e}^{i\theta P}.$$

$${\mathcal{O}}_q=\frac{1}{2L}\sum _{m=0}^{2L-1}{e}^{i\frac{qm\pi }{L}}{e}^{-i\frac{m\pi P}{L}}\mathcal{O}{e}^{i\frac{m\pi P}{L}}.$$

$$I_q=2^{-L}\text{Tr}({\mathcal{O}}_q{\mathcal{O}}_{-q})=2^{-L}\text{Tr}({\mathcal{O}}_q\mathcal{O})=\frac{2^{-L}}{2L}\sum _{m=0}^{2L-1}{e}^{i\frac{qm\pi }{L}}\text{Tr}(e^{-i\frac{m\pi P}{L}}\mathcal{O}{e}^{i\frac{m\pi P}{L}}\mathcal{O})=\frac{1}{2L}\sum _{m=0}^{2L-1}{e}^{i\frac{qm\pi }{L}}{S}_m$$

$$S_m(t)=2^{-L}\text{Tr}(e^{-i\frac{m\pi P}{L}}{U}_t\mathcal{O}(0)U_{-t}{e}^{i\frac{m\pi P}{L}}{U}_t\mathcal{O}(0)U_{-t})=2^{-L}\text{Tr}(\mathcal{O}(0)U_{-t}{e}^{-i\frac{m\pi P}{L}}{U}_t\mathcal{O}(0)U_{-t}{e}^{i\frac{m\pi P}{L}}{U}_t)$$

$$[P,\mathcal{O}]$$

$$[P,[P,\mathcal{O}]]$$

$$[P,[P,[P,\mathcal{O}]]]$$

$$⋮$$

$$\underset{2L}{\underbrace{[P,[P,[P,[P,\cdots [P,}}\mathcal{O}],\cdots ]]]]$$

$$\left[\begin{array}{c}[P,\mathcal{O}]\\ [P,[P,[P,\mathcal{O}]]]\\ [P,[P,[P,[P,[P,\mathcal{O}]]]]]\\ ⋮\end{array}\right]$$

$$\left[\begin{array}{c}[P,[P,\mathcal{O}]]\\ [P,[P,[P,[P,\mathcal{O}]]]]\\ [P,[P,[P,[P,[P,[P,\mathcal{O}]]]]]]\\ ⋮\end{array}\right]$$

$$S_{\varphi }(t_1,t_2)=2^{-L}\text{Tr}[e^{-i\varphi {\mathcal{O}}^{\prime }}Z(t_1)e^{i\varphi {\mathcal{O}}^{\prime }}Z(t_2)]=\sum _q{e}^{-iq\varphi }{I}_q(t_1,t_2)$$

$$2^{-L}\text{Tr}[e^{-i\varphi {\mathcal{O}}^{\prime }}\overline{Z}{e}^{i\varphi {\mathcal{O}}^{\prime }}\overline{Z}]=\sum _q{e}^{-iq\varphi }\frac{1}{t_{\mathrm{pre}}^2}{\int }_0^{t_{\mathrm{pre}}}d{t}_1{\int }_0^{t_{\mathrm{pre}}}d{t}_2{I}_q(t_1,t_2)$$

$$C_{\overline{Z}{\mathcal{O}}^{\prime }}=-4L^{-1}⟨[\overline{Z},{\mathcal{O}}^{\prime }{]}^2{⟩}_{\beta =0}=\sum _q\frac{1}{t_{\mathrm{pre}}^2}{\int }_0^{t_{\mathrm{pre}}}d{t}_1{\int }_0^{t_{\mathrm{pre}}}d{t}_2{I}_q(t_1,t_2)$$

$$C_{\tilde{Z}{\mathcal{O}}^{\prime }}=-4L^{-1}⟨[\tilde{Z},{\mathcal{O}}^{\prime }{]}^2{⟩}_{\beta =0}=\sum _q\frac{q^2}{M^2}\sum _j^M\sum _k^M{I}_q(j\Delta t,k\Delta t)$$

Full text

The dynamics of many-body quantum systems can display a multitude of interesting phenomena, ranging from thermalization Kaufman et al. (2016); Kucsko et al. (2018) to many-body localization (MBL) Basko et al. (2006, 2007); Nandkishore and Huse (2015); Schreiber et al. (2015); Smith et al. (2016); Choi et al. (2016); Wei et al. (2018); Lukin et al. , discrete time crystals Khemani et al. (2016); Keyserlingk et al. (2016); Else et al. (2016); Moessner and Sondhi (2017); Sacha and Zakrzewski (2018); Yao et al. (2017); Zhang et al. (2016); Choi et al. (2017); Ho et al. (2017), and dynamical phase transitions Zhang et al. (2017); Yuzbashyan et al. (2006); Heyl et al. (2013); Žunkovič et al. (2018); Fläschner et al. (2016); Jurcevic et al. (2017). Recently, there has been increased interests in systems exhibiting nonergodic dynamics in the absence of any disorder or incommensurate fields, such as quasi-MBL in translationally invariant systems Yao et al. (2016) and disorder free localization Papić et al. (2015); Smith et al. (2017); Michailidis et al. (2018). Another intriguing possibility is prethermalization, where nonintegrable quantum systems may fail to thermalize on practically accessible timescales Berges et al. (2004); Gring et al. (2012); Else et al. (2017a, b); Abanin et al. (2017); Neyenhuis et al. (2017), due to an emergent quasi-local integral of motion.

Here we study thermalization and prethermalization by measuring out-of-time-ordered (OTO) commutators Larkin and Ovchinnikov (1969); Kitaev (2015); Kitaev and Suh (2017); Garttner et al. (2017); Li et al. (2017); Landsman et al. (2018), which are powerful indicators of information scrambling, but are typically difficult to observe experimentally.

We exploit Hamiltonian engineering techniques to investigate the onset of prethermalization in a nuclear spin system in a natural crystal. We can access different regimes by manipulating the relative strengths of the dipolar interactions among spins and the transverse magnetic field. After a quench, we experimentally measure OTO commutators using multiple quantum coherence (MQC) experiments Munowitz and Pines (1975); Gärttner et al. (2017); Wei et al. (2018) for a system initially at an effective infinite temperature. In the low field regime, the system is thermalizing and the commutator keeps increasing in the observed timescale. In the high field regime, instead, an emergent conserved quantity arises due to prethermalization and the OTO commutator involving such prethermal conserved quantity saturates after a short time. We further support the interpretation of our experimental results by constructing the prethermal Hamiltonian perturbatively Else et al. (2017c); Abanin et al. (2017). We numerically observe the divergence of the perturbation series below a certain transverse field threshold, indicating the breakdown of prethermal dynamics and the onset of the thermal regime.

We work with nuclear spins in fluorapatite (FAp) der Lugt and Caspers (1964), an experimental system recently used to show MBL Wei et al. (2018). The 19F spins-1/2 form linear chains in the crystal and are coupled by the magnetic dipolar interaction. A single crystal is placed in a large (7 T) magnetic field at room temperature. In a strong magnetic field the interaction Hamiltonian for the 19F spins is given by the secular dipolar Hamiltonian $H_{\text{Dipz}}=\sum_{j,k>j}J_{jk}\!\left[S_{z}^{j}S_{z}^{k}-\frac{1}{2}(S_{x}^{j}S_{x}^{k}+S_{y}^{j}S_{y}^{k})\right]$, where $J_{jk}=J|j-k|^{-3}$. Here $S_{\alpha}^{j}$ $(\alpha=x,y,z)$ are spin-1/2 operators of the $j$-th 19F spin (see Section I in SM). In the timescales we explore, the system can be approximately treated as an ensemble of identical spin chains Cappellaro et al. (2007); Zhang et al. (2009); Ramanathan et al. (2011), since the interchain coupling is $\sim\!40$ times weaker than the intrachain coupling. The coupling to 31P spins in the lattice is refocused by the applied control, and the spin-lattice relaxation effects are negligible ($T_{1}\approx 0.8$ s). The dynamics of the 19F spins are thus well approximated by a 1D closed quantum system with dipolar couplings. While the corresponding 1D, nearest-neighbor XXZ Hamiltonian is integrable Alcaraz et al. (1987); Sklyanin (1988); Wang et al. (2016), the Hamiltonian we consider can lead to diffusive Sodickson and Waugh (1995); Zhang and Cory (1998) and chaotic behavior Jyoti (2017) in 3D. In the presence of a transverse field, the system is known to show a quantum phase transition Isidori et al. (2011).

In experiments we consider the dynamics under an engineered Floquet Hamiltonian, obtained by modulating $H_{\text{Dipz}}$ with periodic sequences of strong rf pulse Haeberlen and Waugh (1968) that can also introduce quenches and time reversal. To lowest order of Magnus expansion the pulse sequence (see Section II in SM) engineers a dipolar Hamiltonian along the $y$ direction, $H_{\text{dipy}}$ Wei et al. (2018), while an effective static transverse field is introduced by phase shifting all the pulses. The resulting Floquet-Trotter Hamiltonian is equivalent to its lowest order to a transverse field dipolar Hamiltonian $H_{\text{TDip}}\!=\!uH_{\text{Dipy}}\!+\!gZ$ 111Note that the form of this Hamiltonian is reminiscent of a spin-locking dipolar interaction in NMR Slichter (1996):

[TABLE]

where both $u$ and $g$ are under experimental control (for details see Section II in SM) and we set $J=-uJ_{j,j+1}$ being the engineered nearest-neighbor coupling strength. In all experiments we set $u=0.2$ and pulse sequence period $t_{c}=96~{}\mu$s, which corresponds to an effective $Jt_{c}=0.62$, given the natural $J_{j,j+1}=-33$ krad/s neighbor coupling strength in fluorapatite. Whereas for either $g\!=\!0$ or $J\!=\!0$ the magnetizations $Y=\sum_{j}S_{y}^{j}$ and $Z=\sum_{j}S_{z}^{j}$ are exactly conserved, respectively, for finite values their dynamics is quite different, as observed in the experimental two point correlators $\langle Z(t)Z\rangle_{\beta=0}$ and $\langle Y(t)Y\rangle_{\beta=0}$ shown in Fig. 1(a). Indeed, $H_{\text{TDip}}$ can be mapped to a prethermal Hamiltonian only at high field, while at low field we expect thermalization. As shown in Ref. Abanin et al. (2017), a prethermal regime exists for Hamiltonians that can be divided into two parts $H\!=\!H_{0}+\epsilon V$, with $H_{0}$ having integer eigenvalues up to a scaling factor $C$, $e^{i2\pi CH_{0}}\!=\!\mathbb{1}$. For sufficiently small $\epsilon$, $H$ can be approximately transformed to a prethermal Hamiltonian $H_{\mathrm{pre}}$, through a local unitary $R$ Abanin et al. (2017); Else et al. (2017d), i.e. $RHR^{\dagger}=H_{\text{pre}}+\delta H$, where $\delta H$ is exponentially small in $\epsilon$, $R=\mathbb{1}+O(\epsilon)$, and $[H_{0}^{\prime},H_{\text{pre}}]=0$, where $H_{0}^{\prime}$ in the frame rotated by $R$ has the same matrix representation as $H_{0}$ in the original frame, so they are different physical operators. As the prethermal Hamiltonian conserves $H^{\prime}_{0}$, $R^{\dagger}H_{0}^{\prime}R$ is a conserved quantity in the original frame up to an exponentially long time $t_{\text{pre}}$, after which the small correction $\delta H$ thermalizes the system. In the transverse field dipolar model with $g\gg J$, we can identify the dominant part with the field and the perturbation with the dipolar interaction. The prethermal Hamiltonian is then $H_{\mathrm{pre}}=gZ^{\prime}-uH^{\prime}_{\text{Dipz}}/2+O(J/g)$. Then, in the prethermal regime we expect an emergent conserved quantity, $Z_{\text{pre}}$, related to $Z^{\prime}$ by a local unitary transformation $R$, $Z_{\text{pre}}=R^{\dagger}Z^{\prime}R$.

To investigate the presence of this emergent constant of motion beyond the partial information given by local observables [Fig. 1(a)] we experimentally analyze the properties of OTO commutators [Fig. 1(b-c)], defined as $C_{\mathrm{AB}}(t)\equiv\langle[A(t),B][A(t),B]^{\dagger}\rangle_{\beta}$, where $A(t)=U(t)AU(t)^{\dagger}$, with $U(t)=e^{-i\hat{H}t}$ and $\hat{H}$ being the system Hamiltonian. Here $\langle\cdot\rangle_{\beta}=\mathrm{Tr}(e^{-\beta\hat{H}}\,\cdot\,)/\mathrm{Tr}(e^{-\beta\hat{H}})$ denotes the ensemble average at the inverse temperature $\beta$. The OTO commutator contains a term with an unconventional temporal order, the OTO correlator $F(t)\equiv\langle{A}^{\dagger}(t){B}^{\dagger}A(t)B\rangle_{\beta}$, which can provide a more accurate description of operator scrambling than, e.g., Loschmidt echoes Quan et al. (2006); Swingle et al. (2016); Hahn (1950); Rhim et al. (1971); Zhang et al. (1992); Andersen et al. (2003); Gorin et al. (2006); Prosen et al. (2003). We exploit our ability to engineer a time reversal of the Hamiltonian in Eq. (1) to measure the OTO commutator of extensive observables Kukuljan et al. (2017), as we explain in the following.

In room temperature NMR experiments, the initial state for a chain of $L$ spins is described by the density matrix $\rho(0)\!\approx\!(\mathbb{1}\!-\!\epsilon Z)/2^{L}$, with $\epsilon\!\sim\!10^{-5}$. Since the identity operator does not contribute to any measurable signal, we only care about the deviation from it, $\delta\rho=2Z/\sqrt{L}$, which has been normalized such that $\text{Tr}(\delta\rho^{2})/2^{L}=1$. The mixed initial state enables the experimental study of two-point correlators and OTO commutators in a straightforward way. Since $\delta\rho(0)$ is usually the collective spin magnetization pointing in some direction, $\mathcal{O}_{\bf n}=\sum_{j}{\bf n}\cdot{\boldsymbol{S}}_{j}$ and we can measure the collective magnetization around any axis, the typical signal is the two-point correlation, $4\text{Tr}[U(t)\delta\rho(0)U^{{\dagger}}(t)\mathcal{O}_{\bf n}]/(2^{L}L)\equiv\langle\mathcal{O}_{\bf n}(t)\mathcal{O}_{\bf n}\rangle_{\beta=0}$. That is, in our experiments, the (deviation of) the density matrix plays the role of an observable for an effective simulated system at infinite temperature. Crucially, however, the “simulated observable” $\delta\rho$ will thermalize at long times under the strong driving, $\langle\delta\rho(t)\rangle=0$: this enables distinguishing the prethermal regimes from the expected (zero) signal at long times due to the eventual thermalization. MQC experiments Munowitz and Pines (1975); Baum et al. (1985); Ramanathan et al. (2003) measure the overlap of the time-evolved density matrix, $\delta\rho(t)=U(t)\delta\rho(0)U^{{\dagger}}(t)$, with itself after a collective rotation. The overall measured signal can be expressed as

[TABLE]

Taking a discrete Fourier transform of $S_{\phi}$ with respect to $\phi$ yields the MQC intensities: $S_{\phi}=\sum_{q}e^{-iq\phi}I_{q}$. Expanding $S_{\phi}$ in powers of $\phi$, it can be shown that $\text{Tr}([\delta\rho(t),\mathcal{O}_{\bf n}]^{2})/2^{L}\!=\!-\sum_{q}q^{2}I_{q}$. Setting $\delta\rho(0)=\mathcal{O}_{\bf n^{\prime}}$, we can write

[TABLE]

Eq. (3) is the central idea of our experiments: by measuring the second moment of the MQC intensities encoded in $\delta\rho(t)$ along $\mathcal{O}_{\bf n}$ one can obtain the OTO commutator between $\mathcal{O}_{\bf n^{\prime}}(t)$ and $\mathcal{O}_{\bf n}$ as if the system were at infinite temperature 222Notice that exchanging $\mathcal{O}_{\bf n^{\prime}}$ and $\mathcal{O}_{\bf n}$ will result in a different MQC distribution $I_{q}$, however its second moment remains the same.. Eq. (3) was first derived in a different context in Ref. Khitrin (1997) for NMR systems. When applied to pure states, it relates the second moment of the MQC distribution to the quantum Fisher information Gärttner et al. (2017).

To study the system dynamics after a quench to Hamiltonian (1), we measure the OTO commutator $C_{\mathrm{YZ}}\equiv 4L^{-1}\langle|[Y(t),Z]|^{2}\rangle_{\beta=0}$ for various transverse field strengths and times [see Fig. 1(b-c)]. First note that in the limit $g\rightarrow\infty$, $Z$ is a conserved quantity thus making $C_{\mathrm{YZ}}$ constant. In Fig. 1(c) we observe that for large but finite transverse field $C_{\mathrm{YZ}}$ stops growing at an early time, revealing that $Z$ is approaching the emergent conserved quantity $Z_{\text{pre}}$ - as also indicated by the slow decay and persistent oscillation of the two point correlator $\langle Z(t)Z\rangle_{\beta=0}$ [Fig. 1(a)]. For small transverse field, instead, $C_{\mathrm{YZ}}$ keeps increasing, suggesting that the system is thermalizing [Fig. 1(c)]. We note that in the limit of exactly no transverse field, $Y$ is a conserved quantity thus making $C_{\mathrm{YZ}}$ constant. However, as long as a small field is introduced the system becomes thermal until the field strength induces a transition to a prethermal state: we thus observe a maximum of $C_{\mathrm{YZ}}$ at around $g/J\approx 0.5$ [Fig. 1(b)]. The thermal dynamics for an initial effective infinite temperature state is further indicated by the decay of $\langle Y(t)Y\rangle_{\beta=0}$ in Fig. 1(a) and additional OTO commutators presented below.

To gain further insight into the differences between the thermalizing and prethermal regimes, we measure $C_{\mathrm{ZZ}}$ and $C_{\mathrm{YY}}$, as shown in Fig. 2(a). Because these OTO commutators fluctuate significantly in time, we average them at six different times. As $g$ increases, $Z(t)$ approaches the prethermal conserved quantity $Z_{\text{pre}}$, which itself gets close to $Z$, and $C_{\mathrm{ZZ}}$ gets smaller.

This behavior is only observed for OTO commutators involving at least one operator that overlaps with the emergent conserved quantity, while other commutators, such as $C_{\mathrm{YY}}$, keep growing as if the system were thermal, regardless of the transverse field strength (with the exception of exactly zero field, $g=0$).

While we cannot directly measure $Z_{\mathrm{pre}}$, the time-averaged operator $\overline{Z}=t^{-1}_{\text{pre}}\int_{0}^{t_{\text{pre}}}Z(t)dt$ (where $t_{\text{pre}}$ is the timescale over which the prethermal conserved quantity is present) captures its essential features 333In the limit $t_{\text{pre}}\!\!\rightarrow\!\infty$, $\overline{Z}$ becomes an exact integral of motion and it has also recently been proposed as an experimentally accessible observable for detecting MBL to ergodic phase transition Chandran et al. (2015). Indeed, we can generally write $Z(t)=Z_{\text{pre}}+U(t)(Z-Z_{\mathrm{pre}})U(t)^{\dagger}$: then, in the prethermal regime, the second term is small and fluctuates, yielding $\overline{Z}\approx Z_{\text{pre}}$ after time average. We can approximate $\overline{Z}$ with a discrete time average, $\widetilde{Z}=\sum_{n=1}^{N}Z(t_{n})/N$, by independently varying the forward and backward evolution times in the MQC protocol (see SM for details on the experiments and for a comparison between $\overline{Z}$ and $\widetilde{Z}$). Figure 2(b) shows that $\mathrm{Tr}({\widetilde{Z}}^{2})/\mathrm{Tr}({Z(0)}^{2})\to 1$ as $g$ increases, because the time-varying part of $Z(t)$ is very small for large $g$. Furthermore, $4/L\langle|[\widetilde{Z},Z(0)]|^{2}\rangle_{\beta=0}$ approaches zero at large $g$, suggesting that $\lim_{g\to\infty}\widetilde{Z}=Z$.

To support our interpretation of the experimental results, we explicitly construct the prethermal Hamiltonian, showing that indeed $Z_{\text{pre}}\approx Z$ is an emergent constant of motion.

The prethermal Hamiltonian can be expanded in powers of $\epsilon=J/g$

[TABLE]

and numerically evaluated up to order $n_{M}$ (see Section IV in SM). It has been shown Abanin et al. (2017) that for generic many-body systems the series in Eq. (18) might not converge as $n_{M}\rightarrow\infty$, but there exists an optimal order $n^{*}$ when truncating the series, so that $H_{\text{pre}}$ is most similar to $H$. If the system Hamiltonian does indeed support a prethermal phase, we expect its eigenvalues $E_{m}$ to be close to the prethermal Hamiltonian ones, $E_{m}^{\mathrm{pre}}$. We thus calculate the eigenvalue difference $r\equiv\mathrm{mean}_{m}(E_{m}-E_{m}^{\mathrm{pre}})/L$ (where $m$ labels the eigenvalues in ascending order), expecting $r$ to converge to zero only in the prethermal phase. Figure 3(a) shows $r$ as a function of maximum truncation order $n_{M}$ for different values of $\epsilon$. For large $g$, $r\approx 0$ appears to converge up to the largest numerically accessible order, suggesting that $H_{\mathrm{pre}}$ is similar to $H$ and there exists an approximately conserved quantity $Z_{\mathrm{pre}}$. For small $g$ however, $r$ diverges, indicating that a prethermal Hamiltonian that conserves $Z^{\prime}$ cannot be found. The transition happens at around $g/J=0.5$.

To further demonstrate that a conserved quantity emerges for large $g$, we simulate $Z$ at large times $(Jt=10^{3})$ and decompose it according to the Hamming weight Wei et al. (2018)

[TABLE]

where $\mathcal{B}_{k}^{s}$ are operators composed of tensor products of $k$ Pauli matrices and $L-k$ identity operators, and $\zeta_{k}\propto 3^{k}\times{L\choose k}$ labels the number of configurations with $k$ non-identity Pauli operators. We define the Hamming weight of $k$-spin correlations as $f_{k}=\sum_{s=1}^{\zeta_{k}}[b_{k}^{s}]^{2}$, satisfying $\sum_{k=1}^{L}f_{k}=1$. Figure 3(b) shows that for small transverse field $f_{k}$ is approximately proportional to $\zeta_{k}$, suggesting that all possible operators $\mathcal{B}_{k}^{s}$ have the same weight, in agreement with the eigenvalue thermalization hypothesis Hosur and Qi (2016); Rigol et al. (2008); Deutsch (1991); Srednicki (1994). The result is qualitatively different for $g\gg 1$, where a significant one-body term, $f_{1}$, exists even at very late times ($Jt=10^{3}$), signifying the failure of thermalization and the emergence of the conserved quantity $Z_{\text{pre}}$. We thus study $f_{1}$ as a function of time. For small fields, $g/J\leq 0.5$, the contribution of $f_{1}$ in $Z(t)$ quickly relaxes from one to zero, as shown in Fig. 3(c). For large $g/J$, instead, $f_{1}$ reaches a non-zero, quasi-equilibrium value, signaling the prethermal phase. We do not see the final thermalizing stage in the numerics, possibly because small systems, $L=13$, do not fully thermalize Machado et al. (2017). On the other hand, while $Y$ is conserved at exactly zero field ($g=0$), as soon as a small transverse field is introduced the contribution of $f_{1}$ to $Y(t)$ decays to zero [Fig. 3(d)]. This indicates that the slow dynamics observed for $Y(t)$ at small $g$ is not protected by any prethermal conserved quantity, and will thus thermalize on timescales much shorter than $t_{\text{pre}}$. The quantitative difference between $f_{1}$ for $Y(t)$ and $Z(t)$ can be approximately observed by measuring the two-point correlations $\langle Z(t)Z\rangle_{\beta=0}$ and $\langle Y(t)Y\rangle_{\beta=0}$. As shown in Fig. 1(a), in the small field regime $\langle Y(t)Y\rangle$ decreases rapidly as a function of time, suggesting that $f_{1}(t)$ of $Y(t)$ is not a (prethermal) conserved quantity 444We further numerically demonstrate in the SM that $H$ is unlikely to be in a prethermal regime with respect to NN Ising interactions $\sum_{j}S_{y}^{j}S_{y}^{j+1}$, the dominant term in $H_{\text{TDip}}$ for small values of $g$.. In stark contrast, $\langle Z(t)Z\rangle$ shows a slow decay with periodic oscillations, suggesting that $f_{1}(t)$ of $Z(t)$ is mostly conserved, consistent with prethermalization at large $g/J$ 555taking the average over time, also shows that $\overline{Z}\approx Z$ in the prethermal regime.

In conclusion, we studied the out-of-equilibrium dynamics of the transverse field dipolar interaction in a solid-state NMR quantum simulator. Using MQC techniques, we measured OTO commutators to reveal a distinct dynamics in the high and low field regimes, and identified them as prethermal and thermal phases. In the prethermal regime, when one of the OTO operators is close to the emergent quasi-conserved quantity, the OTO commutator saturates at an early time, while it keeps increasing in the thermal regime, with a transition at about $g/J=0.5$ 666we note that this value does not corresponds to the critical point of the equilibrium quantum phase transition Isidori et al. (2011). We further validate our experimental results numerically, by constructing the prethermal Hamiltonian and verifying the emergence of a conserved quantity at high field, while in the low field regime the dynamics is consistent with thermalization. We demonstrate the value of OTO commutators in investigating non-equilibrium quantum thermodynamics, while also providing a method to experimentally measure OTO commutators that could be extended to other experimental platforms. Similar techniques could be used for example to explore other many-body phenomena, such as localization, dynamics phase transition and information scrambling, paving the way to more comprehensive understanding of out-of-equilibrium quantum many-body systems.

We thank N. Halpern, D. Huse, and I. Cirac for insightful discussions. This work was supported by the National Science Foundation PHY1734011.

SUPPLEMENTARY MATERIAL

I I. Experimental System

The system used in the experiment was a single crystal of fluorapatite (FAp). Fluorapatite is a hexagonal mineral with space group $P6_{3}/m$, with the 19F spin-1/2 nuclei forming linear chains along the $c$-axis. Each fluorine spin in the chain is surrounded by three 31P spin-1/2 nuclei. We used a natural crystal, from which we cut a sample of approximate dimensions 3 mm$\times$3 mm$\times$2 mm. The sample is placed at room temperature inside an NMR superconducting magnet producing a uniform $B=7$ T field. The total Hamiltonian of the system is given by

[TABLE]

The first two terms represent the Zeeman interactions of the F($S$) and P($s$) spins, respectively, with frequencies $\omega_{F}=\gamma_{F}B\approx(2\pi)282.37$ MHz and $\omega_{P}=\gamma_{P}B=(2\pi)121.51$ MHz, where $\gamma_{F/P}$ are the gyromagnetic ratios. The other three terms represent the natural magnetic dipole-dipole interaction among the spins, given generally by

[TABLE]

where $\vec{r}_{ij}$ is the vector between the $ij$ spin pair. Because of the much larger Zeeman interaction, we can truncate the dipolar Hamiltonian to its energy-conserving part (secular Hamiltonian). We then obtain the homonuclear Hamiltonians

[TABLE]

and the heteronuclear interaction between the $F$ and $P$ spins,

[TABLE]

with $J_{jk}=\hbar\gamma_{j}\gamma_{k}\frac{1-3\cos(\theta_{jk})^{2}}{|\vec{r}_{jk}|^{3}}$, where $\theta_{jk}$ is the angle between the vector $\vec{r}_{jk}$ and the magnetic field $z$-axis. The maximum values of the couplings (for the closest spins) are given respectively by $J^{F}=-32.76$ krad s*-1*, $J^{P}=1.20$ krad s*-1* and $J^{FP}=6.12$ krad s*-1*.

The dynamics of this complex many-body system can be mapped to a much simpler, quasi-1D system. First, we note that when the crystal is oriented with its $c$-axis parallel to the external magnetic field the coupling of fluorine spins to the closest off-chain fluorine spin is $\approx 40$ times weaker, while in-chain, next-nearest neighbor couplings are $8$ times weaker. Previous studies on these crystals have indeed observed dynamics consistent with spin chain models, and the system has been proposed as solid-state realizations of quantum wires  Cappellaro et al. (2007, 2011); Ramanathan et al. (2011). This approximation of the experimental system to a 1D, short-range system, although not perfect has been shown to reliably describe experiments for relevant time-scales Rufeil-Fiori et al. (2009); Zhang et al. (2009). The approximation breaks down at longer time, with a convergence of various effects: long-range in-chain and cross chain couplings, as well as pulse errors in the sequences used for Hamiltonian engineering. In addition, the system also undergoes spin relaxation, although on a much longer time-scale ($T_{1}=0.8~{}$s for our sample).

I.1 Error analysis

The desired quantity $S_{m}=\mathrm{Tr}(\delta\rho(t_{\mathrm{end}})\mathcal{O})$, where $\delta\rho(t_{\mathrm{end}})=U_{-t}e^{-i\frac{m\pi P}{L}}U_{t}\delta\rho(0)U_{-t}e^{i\frac{m\pi P}{L}}U_{t}$ is the nontrivial part of the density matrix at the end of the entire pulse sequence $t_{\mathrm{end}}$. To get the uncertainty of $S_{m}$, instead of just measuring one point, we continuously monitor the free evolution of $\delta\rho(t_{\mathrm{end}})$ under the natural Hamiltonian $H_{\mathrm{dipz}}$, so called free induction decay (FID, a typical FID process is shown in Fig. 5). We take the standard deviation of the last 20 data points in the FID as the uncertainty of the $S_{m}$, and then linearly propagate to get the error bar of the OTO commutators.

II II. Hamiltonian Engineering

The Hamiltonians used in the main text were obtained stroboscopically (Floquet Hamiltonians) by applying periodic rf pulse trains to the natural dipolar Hamiltonian that describes the system. We used Average Hamiltonian Theory (AHT Haeberlen and Waugh (1968)) as the basis for our Hamiltonian engineering method, to design the control sequences and determine the approximation errors.

To see how repeatedly applying a periodic pulse sequence modifies the dynamics of the system, we write the total Hamiltonian as $H=H_{\text{dip}}+H_{\text{rf}}$, where $H_{\text{dip}}=\frac{1}{2}\sum_{j<k}J_{jk}(2S_{z}^{j}S_{z}^{k}-S_{x}^{j}S_{x}^{k}-S_{y}^{j}S_{y}^{k})+\sum_{j}h_{j}S_{z}^{j}$ is the system Hamiltonian, and $H_{\text{rf}}(t)$ is the external Hamiltonian due to the rf-pulses. The density matrix $\rho$ evolves under the total Hamiltonian according to $\dot{\rho}=-i[H,\rho]$. Consider an interaction frame defined by $\rho^{\prime}={U_{\text{rf}}}^{\dagger}\rho U_{\text{rf}}$, where $U_{\text{rf}}(t)=\mathcal{T}\exp[-i\int_{0}^{t}H_{\text{rf}}(t^{\prime})dt^{\prime}]$ and $\mathcal{T}$ is the time ordering operator. In this toggling frame, $\rho^{\prime}$ evolves according to $\dot{\rho}^{\prime}=-i[H(t),\rho^{\prime}]$, where $H(t)={U_{\text{rf}}}^{\dagger}H_{\text{dip}}U_{\text{rf}}$. Since $U_{\text{rf}}$ is periodic, $H(t)$ is also periodic with the same period $t_{c}$. The evolution operator over one period can be written as $U(t_{c})=\exp[-iH_{\text{Flq}}t_{c}]$, where $H_{\text{Flq}}$ is called the Floquet Hamiltonian (or in the language of NMR the Average Hamiltonian). Note that if the pulse sequence satisfies the condition $U_{\text{rf}}(t_{c})=1$, the dynamics of $\rho$ and $\rho^{\prime}$ are identical when the system is viewed stroboscopically, i.e., at integer multiples of $t_{c}$. The system evolves as if under a time-independent Hamiltonian $H_{\text{Flq}}$. To calculate $H_{\text{Flq}}$ we employ the Magnus expansion as is usual in AHT: $H_{\text{Flq}}=H^{(0)}+H^{(1)}+\cdots$. The first two terms are given by

[TABLE]

The zeroth order of the average Hamiltonian $H^{(0)}$ is often a good approximation to the Floquet Hamiltonian $H_{\text{Flq}}$, as the first order can be set to zero by simple symmetrization of the pulse sequence.

The basic building block of the pulse sequences we used in this work is given by a 4-pulse sequence Kaur and Cappellaro (2012); Yen and Pines (1983) originally developed to study MQC. We denote a generic 4-pulse sequence as $P(\tau_{1},{\bf n}_{1},\tau_{2},{\bf n}_{2},\tau_{3},{\bf n}_{3},\tau_{4},{\bf n}_{4},\tau_{5})$, where ${\bf n}_{j}$ represents the direction of the $j$-th $\pi/2$ pulse, and $\tau_{j}$’s the delays interleaving the pulses. In our experiments, the $\pi/2$ pulses have a width $t_{w}$ of typically 1 $\mu$s. $\tau_{j}$ starts and/or ends at the midpoints of the pulses (see also Fig. 4). In this notation, our forward 16-pulse sequence can be expressed as

[TABLE]

and the backward sequence as

[TABLE]

where $\{{\bf\bar{x}},{\bf\bar{y}}\}\equiv\{{\bf-x},{\bf-y}\}$. The delays are given by

[TABLE]

where $\tau$ is typically 4 $\mu$s. The cycle time $t_{c}$, defined as the total time of the sequence, is given by $t_{c}=24\tau$. $u$ is a dimensionless adjustable parameter, and is restricted such that none of the inter-pulse spacings becomes negative. We set $u=0.2$ in experiments.

Our forward and backward pulse sequences with finite pulse width produce average dipolar Hamiltonians with opposite sign:

[TABLE]

and $H^{(1)}=0$ (the first order can always be set to zero by a proper symmetrization of the sequence). The interaction between $F$ and $P$ spins averages to zero under this pulse sequence. The average Hamiltonian corresponding to the forward sequence is presented as Eq.(1) in the main text.

A uniform transverse field can be introduced in $H^{(0)}$ by phase-shifting the entire pulse sequence. Consider rotating the $n$-th cycle of the pulse sequence by $(n-1)\phi$ around the ${\bf z}$ axis, which can be accomplished by phase shifting all the pulse directions ${\bf n}_{j}$ in the $n$-th cycle by $(n-1)\phi$. The evolution operator for each cycle is given by

[TABLE]

where $Z=\sum_{j}S_{z}^{j}$. The total evolution operator over $n$ cycles is given by the product:

[TABLE]

where the total Hamiltonian is given by $H_{\text{t}}=H_{\mathrm{Flq}}+gZ$, with $g=-\phi/t_{c}$. We note that for the values we considered, the approximation $\left[e^{i\phi Z}e^{-iH_{\mathrm{Flq}}t_{c}}\right]^{n}\approx e^{-i\left(H^{(0)}-\frac{\phi}{t_{c}}Z\right)nt_{c}}$ is extremely good, with the error being on order $g/J\sim 1.5\times 10^{-3}$. The rotation approach also generates an extra term $e^{-in\phi Z}$, this term can be canceled in MQC experiments by rotating the encoding pulse by $n\phi$. In the backward sequence, the phase shifting is done in reverse order. We note that the phase shift is typically not periodic over the course of the experiment (for example, for small values of $g$, even the maximum $n\phi<2\pi$).

We note that our methods can be applied more broadly to engineer desired Hamiltonians $H_{des}$ using only collective rotations of the spins applied to the naturally occurring Hamiltonian, $H_{nat}$. The engineered Hamiltonian is obtained by piece-wise constant evolution under rotated versions of the natural Hamiltonian under the condition $\sum_{k}R_{k}H_{nat}R_{k}^{\dagger}=H_{des},$ where $R_{k}$ are collective rotations of all the spins, which achieves the desired operator to first order in a Magnus expansion. Symmetrization of the sequence can further cancel out the lowest order correction. Using only collective pulses limits which Hamiltonians can be engineered, due to symmetries of the natural Hamiltonian and the action of collective operators. For typical two-body interactions of spin-1/2, an efficient tool to predict which Hamiltonian are accessible is to use spherical tensors Ajoy and Cappellaro (2013).

III III. Multiple Quantum Coherence

Multiple quantum coherence (MQC) can be understood as a Fourier decomposition of a many-spin operator $\mathcal{O}$ with respect to another operator $\mathcal{P}$ satisfying $e^{-2im\pi P}=\mathbb{1}$, where $m$ is an integer. For simplicity consider the case where $P=\sum_{j}{\bf S}_{j}\cdot{\bf n}$ is the generator of global spin rotation around the ${\bf n}$-axis. Any traceless hermitian operator $\mathcal{O}$ can be decomposed as

[TABLE]

where $L$ is the total number of spins in the system, and $\mathcal{O}_{q}$ satisfies following two equivalent relations:

[TABLE]

From the first relation it is evident that $\mathcal{O}_{-q}=\mathcal{O}_{q}^{\dagger}$. For a given $\mathcal{O}$, each multiple quantum component $\mathcal{O}_{q}$ can be formally written as

[TABLE]

which indicates that the multiple quantum components $\mathcal{O}_{q}$ come from the Fourier decomposition of $\mathcal{O}$ with respect to $P$. In addition, it can also be shown that different $\mathcal{O}_{q}$ are orthogonal: $\text{Tr}(\mathcal{O}_{q}\mathcal{O}_{q^{\prime}})\propto\delta_{q(-q^{\prime})}$. Using the fact that the maximum $q$ allowed is equal to the total number of spins $L$, one can replace Eq. (12) as a discrete Fourier transform:

[TABLE]

The MQC intensities are defined as $I_{q}=2^{-L}\text{Tr}(\mathcal{O}_{q}\mathcal{O}_{-q})$. Using orthogonality condition and Eq. (13) it can be shown that

[TABLE]

where $S_{m}=2^{-L}\text{Tr}(e^{-i\frac{m\pi P}{L}}\mathcal{O}e^{i\frac{m\pi P}{L}}\mathcal{O})$ is the signal we measure in the experiments (Eq. (2) in the main text). Measuring $S_{m}$ requires both forward and backward time evolution. Since $\mathcal{O}(t)=U_{t}\mathcal{O}(0)U_{-t}$, the signal at time $t$ can be written as

[TABLE]

where we used the cyclic property of the trace in the second step. The MQC protocol consists of four steps: forward time evolution by $U_{t}$, rotation (phase tagging) by $e^{i\frac{m\pi P}{L}}$, backward time evolution by $U_{-t}$, and detection of initial state $\mathcal{O}(0)$. A schematic of the experimental protocol is illustrated in Fig. 4.

Eq. (14) gives us a way to find $I_{q}$ experimentally without knowing the individual $\mathcal{O}_{q}$. One can use the second relationship in Eq. (10) and (11) to find $\mathcal{O}_{q}$ explicitly. Starting from $[P,\mathcal{O}]=\sum_{q}q\mathcal{O}_{q}$, by repeatedly commuting both sides with $P$ one can generate a linear system in $\mathcal{O}_{q}$:

[TABLE]

Using $\mathcal{O}_{q}^{\pm}=\mathcal{O}_{q}\pm\mathcal{O}_{-q}$, one can rewrite the above into two linear systems of equations:

[TABLE]

The constant matrices are in the form of Vandermondes matrix, and can be exactly inverted. Thus, by calculating nested commutators of $\mathcal{O}$ and $P$ one can find all the multiple quantum components $\mathcal{O}_{q}$. In practice we find Eq. (LABEL:eq:mqcfind) to be more efficient at finding $\mathcal{O}_{q}$ than Eq. (13).

III.1 OTO commutators involving the prethermal integrals of motion

We are interested in measuring the OTO commutator $C_{\overline{Z}\mathcal{O^{\prime}}}=-4L^{-1}\langle[\overline{Z},\mathcal{O^{\prime}}]^{2}\rangle_{\beta=0}$, where the prethermal integral of motion $\overline{Z}$ is given by the time average of operator $Z(t)$: $\overline{Z}=\frac{1}{t_{\mathrm{pre}}}\int_{0}^{t_{\mathrm{pre}}}Z(t)$, where $t_{\mathrm{pre}}$ is the timescale over which the prethermal conserved quantity is present. This commutator can be found by measuring the MQC intensities with different forward and backward evolution times. The measured signal in MQC experiments is given by

[TABLE]

where $t_{1}$ and $t_{2}$ are now independent. Integrate both sides with respect to $t_{1}$ and $t_{2}$ upto $t_{\mathrm{pre}}$ gives

[TABLE]

expand both sides with respect to $\phi$ and equate terms in $\phi^{2}$ gives

[TABLE]

In the experiments we can only measure $I_{q}(t_{1},t_{2})$ at discrete time points. If we replace the integrals by discrete sums we can measure

[TABLE]

where $\tilde{Z}=\frac{1}{M}\sum_{j}Z(j\Delta t)$ is a discrete time average of $Z(t)$. The prethermal time scale $t_{\mathrm{pre}}$ is very large, so we assume the experimentally measured time is smaller than $t_{\mathrm{pre}}$. $\tilde{Z}$ can be used as an approximation for the integral of motion $\overline{Z}$.

IV IV. Constructing the Prethermal Hamiltonian

In this section we review the mathematical framework for prethermalization in time-independent systems. The derivations are based on Abanin et al. Abanin et al. (2017), but using a different approach to calculate the higher order correlations to the prethermal Hamiltonian. In prethermal systems, certain observables can exhibit slow, non-ergodic, dynamics even in the absence of any disorder. The reason is that the prethermal Hamiltonian has an emergent symmetry given by $[H_{\text{pre}},H_{0}^{\prime}]=0$ which is preserved for an exponentially long prethermal time. Thus, any operator that overlaps with $H_{0}^{\prime}$ will have slow dynamics due to the emergent symmetry. The exponential timescale for $t_{\text{pre}}$ has been proved rigorously for both time-dependent and independent systems in Abanin et al. (2017).

To begin, it is useful to divide the transverse field dipolar Hamiltonian into two parts: $H_{\text{TDip}}=H_{0}+V$ up to a scaling factor, where

[TABLE]

where $\epsilon$ is a small parameter. The idea is to perform a local unitary transformation to bring the Hamiltonian into the following form

[TABLE]

where $[D,H^{\prime}_{0}]=0$, $[\delta H,H_{\text{pre}}]\neq 0$, and $H_{\text{pre}}=H_{0}^{\prime}+D$. If $\delta H$ is zero, then $R^{\dagger}H_{\text{pre}}R$ is an exact conserved quantity. Any operators not orthogonal to $R^{\dagger}H_{\text{pre}}R$ will have some finite component conserved up to infinite times, and the prethermal time $t_{\text{pre}}$ is infinity. If, however, $\delta H$ is nonzero but exponentially small compared to $H_{\text{pre}}$, then $R^{\dagger}H_{\text{pre}}R$ is an emergent symmetry conserved only up to $t_{\text{pre}}$.

In order to calculate the prethermal effective Hamiltonian $H_{\text{pre}}$, let $R=e^{S}$ and expand $S$, and $D$ in powers of $\epsilon$:

[TABLE]

Once $D_{j}$ are found the prethermal effective Hamiltonian up to $n_{M}$-th order is given by

[TABLE]

It is not known whether the series in Eq. (18) will converge for any $\epsilon$ when $n_{M}\rightarrow\infty$. When the series does converge then prethermal Hamiltonian is an exact symmetry.

It is instructive to see how to solve this perturbation series iteratively. Using the Baker-Campbell-Hausdorff formula

[TABLE]

on the LHS of Eq. (17) and collecting terms with same order of $\epsilon$, it can be shown that the equations up to $\epsilon^{5}$ are

[TABLE]

where $i\leftrightarrow j$ indicates permutation of $S_{i}$ and $S_{k}$ in the commutators. Starting at first order in $\epsilon$, all terms on the RHS besides $h_{j}\equiv[S_{j},H_{0}]$ are split into two parts: the part that commutes with $H_{0}$ is stored in $D_{j}$; $S_{j}$ is then chosen such that $[S_{j},H_{0}]$ exactly cancels the other part, which does not commute with $H_{0}$, thus automatically making $V^{\prime}_{j}=0$. The iteration is repeated in the next order.

Since $H_{0}$ satisfies $e^{2im\pi H_{0}}=\mathbb{1}$, one can use the same approach for calculating MQC components to systematically and uniquely find $D_{j}$ and $S_{j}$. At each order, decompose $h_{j}$ to its MQC components: $h_{j}=\sum_{q}h_{jq}$, where $[H_{0},h_{jq}]=qh_{jq}$, and find them using Eq. (LABEL:eq:mqcfind). In terms of $h_{jq}$, $D_{j}$ and $S_{j}$ can be written as

[TABLE]

The MQC components $h_{jq}$ are referred to as generalized ladder operators in Lin and Motrunich (2017). If the original Hamiltonian contains only 2-body terms, $S_{j}$ contains at most $(j+1)$-body operator, guaranteeing the locality of the transformation.

To confirm the locality numerically, we decompose $S$ to sum of products of Pauli operators. Each product operator is assigned a correlation distance defined as the distance between the furtherest non-identity operators (e.g. both $\sigma_{0}\otimes\sigma_{1}\otimes\sigma_{0}\otimes\sigma_{2}\otimes\sigma_{0}$ and $\sigma_{0}\otimes\sigma_{1}\otimes\sigma_{3}\otimes\sigma_{2}\otimes\sigma_{0}$ are of correlation distance 3). The weight of terms of different correlation distance is shown in Fig. 6. The local-to-nonlocal transition happens at around $g/J=0.5$. The problem of finding prethermal Hamiltonian is intimately related to finding the slowest operators, or the operators that best commutes with the Hamiltonian Kim et al. (2015); Lin and Motrunich (2017).

IV.1 No Prethermalization for small $g$

In Fig. 7 we explore whether the Hamiltonian in Eq. (1) can be prethermal for small values of $g$. In the main text, it is shown that a prethermal Hamiltonian cannot be generated using $H_{0}=Z$ for small $g$. Another possible prethermal generator is $H_{0}=\sum_{j}2S_{y}^{j}S_{y}^{j+1}$. Following the method given in the previous section, we construct the prethermal Hamilontian order by order, and calculate the mean eigenvalue difference between the prethermal Hamiltonian and the original Hamiltonian. The eigenvalue difference $r$ stops decreasing after $n_{M}>2$ even for moderate transverse field strength (Fig. 7), while for large transverse field it diverges as expected. The prethermalization theory fails to apply to small field regime because even with only dipolar interaction $||H_{0}||\approx||V||$, thus the series expansion does not converge.

V V. Extended Data

V.1 Experimental data for $C_{YY}(t)$ and $C_{ZZ}(t)$

Figure 8 shows the experimentally measured $C_{YY}(t)$ and $C_{ZZ}(t)$ at six different times, whose averaged value is presented in Fig. 1(c) in the main text.

V.2 Experimental data for $\text{Tr}(Z(t)Z)$ and $\text{Tr}(Y(t)Y)$

The extended data for $\text{Tr}(Z(t)Z)$ and $\text{Tr}(Y(t)Y)$ is shown in Fig. 9. The $\text{Tr}(Z(t)Z)$ exhibits oscillations on top of a decay. After averaging over one oscillation period, $\text{Tr}(Z(t)Z)$ can be fitted to an exponential function $Ae^{-\gamma Jt}$, where $\gamma$ is the dimensionless decay rate. Figure 9(c) shows the decay rate $\gamma$ of $\text{Tr}(Z(t)Z)$ can be approximated by an exponential function of field strength $J/g$, i.e. $\gamma=\gamma_{0}e^{-\alpha g/J}+\gamma_{\infty}$, where $\gamma_{\infty}$ is a background decay in the experiments that is independent of $g$. This exponentially slow decay agrees with the theory for Floquet prethermalization Else et al. (2017a, b); Abanin et al. (2017). $\text{Tr}(Y(t)Y)$ shows exponential decay at small $g/J$ [Fig. 9(b)], whose decay rate is shown in Fig. 9(d). The scaling for the decay rate of $\text{Tr}(Y(t)Y)$ cannot be extracted due to limited data points.

VI VI. Numerical results

In this section we show the numerical results of OTO commutators using exact diagonalization (Fig. 10). System size $L=12$ and open boundary condition are used here. In comparison with Fig. 1 and Fig. 2 in the main text (reproduced here in (e)-(h) in Fig. 10), the numerical and experimental results show quantitatively similar behavior. The difference is mainly due to experimental imperfections, decoherence and finite MQC encodings, as well as the finite size effect in simulation . In Fig. 10(d), we calculated $Z_{\text{inf}}$ exactly by taking the diagonal ensemble in the eigenbasis: $Z_{\text{inf}}=\sum_{n}Z_{nn}\outerproduct{n}{n}$, where $\ket{n}$ is an eigenvector of the Hamiltonian. The measured data in Fig. 10(h) overestimates Fig. 10(d) due to finite discrete averaging.