Detection of quantum phases via out-of-time-order correlators

TL;DR

This paper explores how out-of-time-order correlators (OTOCs) can be used to detect quantum phase transitions, providing an analytical framework and numerical validation with the XXZ model.

Contribution

It introduces a new analytical method linking OTOC dynamics to quantum phases, applicable to various quantum phase transitions.

Findings

OTOCs can dynamically indicate quantum phases.

Analytical relations between OTOCs and phase transitions are established.

Numerical confirmation with the XXZ model supports the theory.

Abstract

We elucidate the relation between out-of-time-order correlators (OTOCs) and quantum phase transitions via analytically studying the OTOC dynamics in a degenerate spectrum. Our method points to key ingredients to dynamically detect quantum phases via out-of-time-order correlators for a wide range of quantum phase transitions and explains the existing numerical results in the literature. We apply our method to a critical model, XXZ model that numerically confirms our predictions.

Equations

$$F(t)$$

$$F(t)=\sum _{\alpha ,\beta ,\gamma ,{\gamma }^{\prime }}{c}_{\alpha }^{\ast }{b}_{\beta }{e}^{-i(E_{\beta }-E_{\alpha }+E_{\gamma }-E_{{\gamma }^{\prime }})t}{W}_{\alpha \gamma }^{†}{V}_{\gamma {\gamma }^{\prime }}^{†}{W}_{{\gamma }^{\prime }\beta },$$

$$F_{t\to \infty }$$

$$F(t\to \infty )$$

$$F_{\text{gs}}(t\to \infty )$$

$$H$$

$$F(t)=\sum _{\alpha ,\beta ,\gamma ,{\gamma }^{\prime }}{c}_{\alpha }^{\ast }{b}_{\beta }{e}^{-i(E_{\beta }-E_{\alpha }+E_{\gamma }-E_{{\gamma }^{\prime }})t}{W}_{\alpha \gamma }^{†}{V}_{\gamma {\gamma }^{\prime }}^{†}{W}_{{\gamma }^{\prime }\beta },$$

$$F_{gs}(t)$$

$$P_{\alpha }$$

$$\text{Tr}(V\mathbb{I})\sim \text{Tr}(V|\psi (0)⟩⟨\psi (0)|)=⟨\psi (0)|V|\psi (0)⟩$$

$$\displaystyle P\left[\bigg{|}\text{Tr}V-\Braket{\psi}{V}{\psi}\bigg{|}\geq\epsilon\right]\lesssim\text{exp}\left(-\mathcal{O}(d_{u})\right),$$

Full text

††preprint: APS/123-QED

Out-of-time-order correlators (OTOCs) Larkin and Ovchinnikov (1969) probe information scrambling in quantum systems of different nature Sekino and Susskind (2008); Lashkari et al. (2013); Swingle and Chowdhury (2017); Hashimoto et al. (2017); Gärttner et al. (2017); Chen et al. (2017); Bohrdt et al. (2017); Daǧ and Duan (2019) and reflect the symmetries Chen et al. (2017); Swingle and Chowdhury (2017); Daǧ and Duan (2019) or lack thereof Sekino and Susskind (2008); Maldacena et al. (2016); Chen et al. (2017) of the underlying Hamiltonian. An OTOC, unlike a time-ordered four-point (or two-point) correlator Huang et al. (2017), can determine the spatial and temporal correlations throughout the system, thus giving rise to a bound on information spread Luitz and Bar Lev (2017); Xu and Swingle (2018); Daǧ and Duan (2019). Through such bounds and the decay rate of an OTOC, one can dynamically detect thermal Chen et al. (2017); Bohrdt et al. (2017); Daǧ and Duan (2019) and localized phases Chen et al. (2017); Swingle and Chowdhury (2017); Huang et al. (2017); Fan et al. (2016); He and Lu (2017); Daǧ and Duan (2019). Recently OTOC has been numerically observed to be susceptible to phase transitions either signaling criticality in a diverging Lyapunov exponent Shen et al. (2017) or showing signatures of symmetry-broken phases in its saturation value Heyl et al. (2018). The latter led to more research that shows the relation emerging in other forms, e.g. in excited states Wang and Pérez-Bernal (2018), or with more experimentally-relevant platforms and system parameters Sun et al. (2018). The interest in providing more verification for such an emergent relation is understandable, not only because the relation points to a practical potential for OTOC in dynamically probing quantum criticality, but also the relation is received as unexpected Heyl et al. (2018). It is indeed an intriguing question how a chaos-detecting and out-of-time ordered correlator that is contributed by presumably all spectrum could also probe ground state physics. The reasons of this relation remain unknown as well as an answer to whether the relation is universal. Motivated by these questions, here, we develop a method on OTOC dynamics to obtain intuition for the emerging relation between quantum phase transitions and out-of-time-order correlators. Remarkably it is possible to dynamically decompose OTOC and show that the ground state physics is the leading order contribution to it under the criteria that our method provides. This is the origin why OTOC saturation value could detect the ground state degeneracy. Therefore, we reach to the conclusion that the OTOC is susceptible to long-range order, while the quasi-long range order is not visible to it. Our method provides additional insights regarding the relation, e.g. (i) the relation is not restricted to already-studied models and 1D Heyl et al. (2018); Sun et al. (2018); (ii) the relation can be extended to include the phase transitions in other eigenstates Wang and Pérez-Bernal (2018). Hence, our theory elucidates the reasons of this unexpected connection, renders it intuitive and universal with further insights. To verify our method, we study the dynamics of 1D critical XXZ chain, where there are Ising and critical XY phases.

Method. Our aim is to be able to come up with an expression that predicts the saturation value of OTOC for long times in the spirit of Eigenstate Thermalization Hypothesis (ETH) Srednicki (1999, 1995). The out-of-time-order correlation function can be defined as

[TABLE]

where $V$ and $W$ are local operators and the expectation value is over an initial state $\Ket{\psi(0)}$. This initial state could be chosen as the ground state Heyl et al. (2018); Gärttner et al. (2017), or a random Haar-distributed state Luitz and Bar Lev (2017); Daǧ and Duan (2019) to approximate an equiprobable state $\mathcal{I}$ in Eq. (1) Sugiura and Shimizu (2012); Luitz and Lev (2017); sup . Eventually, the original definition that is the commutator growth $-\text{Tr}\left(\frac{\text{exp}\left[-\beta H\right]}{Z}\left[W(t),V\right]^{2}\right)$ Maldacena et al. (2016) could be reexpressed in terms of the OTOC of operators $W$ and $V$ with an initial state at the inverse temperature $\beta$. Therefore we can probe the information scrambling through OTOCs Li et al. (2017); Gärttner et al. (2017); Swingle et al. (2016); Bohrdt et al. (2017); Daǧ and Duan (2019).

Given $\Ket{\psi(t)}=\sum_{\alpha}c_{\alpha}e^{-iE_{\alpha}t}\Ket{\psi_{\alpha}}$, where $\Ket{\psi_{\alpha}}$ are eigenstates of the Hamiltonian with the associated eigenvalues $E_{\alpha}$, we define a modified initial state $\Ket{\psi^{\prime}(0)}=V\Ket{\psi(0)}$ and have $\Ket{\psi^{\prime}(t)}=\sum_{\beta}b_{\beta}e^{-iE_{\beta}t}\Ket{\psi_{\beta}}$. Then the OTOC, Eq. (1), can be recast to a fidelity measure of 3-point function, and with the help of completeness relation $\sum_{\gamma}\Ket{\psi_{\gamma}}\Bra{\psi_{\gamma}}=\mathbb{I}$ becomes

[TABLE]

where $\Bra{\psi_{\alpha}}W\Ket{\psi_{\gamma}}=W_{\alpha\gamma}$ are eigenstate expectation values Rigol et al. (2008). Now one can derive the saturation value for long times as well as dynamical features, such as revival timescales in integrable Hamiltonians Dağ et al. (2018).

We study the saturation value in long times, since this value is expected to contain the signature of quantum phases. For long enough times, equilibration in OTOC dynamics can be obtained only when the phase decoheres. Then the equilibration value can be obtained by requesting $E_{\beta}-E_{\alpha}+E_{\gamma}-E_{\gamma^{\prime}}=0$. This condition can be satisfied with four different scenarios: (i) $E_{\alpha}=E_{\beta}$ and $E_{\gamma}=E_{\gamma^{\prime}}$; (ii) $E_{\alpha}=E_{\gamma}$ and $E_{\beta}=E_{\gamma^{\prime}}$; (iii) $E_{\alpha}=E_{\beta}=E_{\gamma}=E_{\gamma^{\prime}}$, which is contained both in (i) and (ii); (iv) $E_{\beta}-E_{\alpha}+E_{\gamma}-E_{\gamma^{\prime}}=0$ with $E_{\beta}\neq E_{\alpha}\neq E_{\gamma}\neq E_{\gamma^{\prime}}$. If a nondegenerate spectrum is assumed, i.e. $E_{\alpha}=E_{\beta}$ implies $\alpha=\beta$, the OTOC reduces to,

[TABLE]

with four terms corresponding to four conditions (i)-(iv), respectively. We note that writing OTOC as in Eq. (Detection of quantum phases via out-of-time-order correlators) proved to be useful previously to understand the quantum chaotic systems better, e.g. in chaotic spin chains with conserved quantities that also obey ETH, decay to zero is not supposed to be exponential, but inverse polynomial in system size Huang et al. (2019) and OTOCs capture eigenstate correlations that ETH cannot Chan et al. (2019). These correlations can readily be seen in the first, second and the fourth terms of Eq. (Detection of quantum phases via out-of-time-order correlators). See Supplement S5 sup for some remarks that immediately follow from Eq. (Detection of quantum phases via out-of-time-order correlators) about systems with nondegenerate chaotic spectra. Now we are going to generalize Eq. (Detection of quantum phases via out-of-time-order correlators) to a more generic form, which allows degeneracy in the energy spectra, because a quantum phase transition usually involves energy degeneracy, e.g. degeneracy from spontaneous symmetry breaking or other sources Dağ et al. (2019). We group all eigenstates of the Hamiltonian into degenerate sets labeled by $\theta$, and each state in its corresponding set is denoted by $\alpha$ for an eigenstate $\psi_{[\theta,\alpha]}$. The OTOC can be reorganized with the new notation, which is one main result of this Letter,

[TABLE]

Here, $\theta,\theta^{\prime},\phi,\phi^{\prime}$ denote degenerate sets while $\alpha,\beta,\gamma,\gamma^{\prime}$ denote quantum states in their corresponding sets. Eq. (Detection of quantum phases via out-of-time-order correlators) can predict the saturation value of OTOC accurately if the OTOC saturates at a finite time. If the OTOC does not saturate or shows transient effects, Eq. (Detection of quantum phases via out-of-time-order correlators) still predicts the time-average of OTOC signal $\bar{F}=1/\mathcal{T}\int dtF(t)$ over a time interval $\mathcal{T}$ with sufficient accuracy. In this sense, Eq. (Detection of quantum phases via out-of-time-order correlators) is not limited to long-time dynamics $t\rightarrow\infty$ sup .

We look for the criteria that the ground state subspace contribution is leading order in the OTOC saturation value Eq. (Detection of quantum phases via out-of-time-order correlators). For this, we first set $W=V$ as the order parameter operator in Eq. (Detection of quantum phases via out-of-time-order correlators) for convenience. Then we expand the coefficients $b_{[\theta,\beta]}=\sum_{\kappa,\tau}W_{[\theta,\beta][\kappa,\tau]}c_{[\kappa,\tau]}$ in Eq. (Detection of quantum phases via out-of-time-order correlators) by using the initial state. If (i) the initial state is set to the state where the phase transition is expected to happen, e.g. the ground state(s) $c_{[1,1]}=1$; and (ii) we apply an ansatz on the matrix elements of the operator projected on this state, e.g. $|W_{[1,\alpha][\theta,\beta]}|^{2}\ll 1$, where $\theta\neq 1$ is a different energy subspace than the subspace of the ground state(s), we observe the following dynamical decomposition:

[TABLE]

Here $F_{\text{gs}}(t\rightarrow\infty)$ is the ground subspace contribution, whereas the $F_{\text{ex}}(t\rightarrow\infty)$ is the contribution of higher energy excitations. The latter is a correction to the ground-state physics in the OTOC, when the criteria are satisfied. The assumption on the initial state sets the scrambling discussed in the rest of the paper to effectively zero temperature. Whereas the operator ansatz becomes even more specific for the phase of interest. If there is a symmetry-broken long-range order to capture, the fluctuations between the matrix elements of the operator are suppressed in the ground state subspace, meaning there is at least a pair of matrix elements accumulating the order $\rightarrow|W_{[1,\alpha][1,\beta]}|^{2}\sim\mathcal{O}\left(1\right)$. This modifies the operator ansatz as $|W_{[1,\alpha][1,\beta]}|^{2}\gg|W_{[1,\gamma][\theta,\gamma^{\prime}]}|^{2}$ for the ordered phase. Thus, we derive the expression for $F_{\text{gs}}(t\rightarrow\infty)$ in the ordered phase as,

[TABLE]

while the operator ansatz simultaneously implies that the OTOC is dominated by the ground state physics, $F_{\text{gs}}\gg F_{\text{ex}}$ in the ordered phase. On the other hand, the fluctuations between the matrix elements of the operator are maximal in a disordered phase, implying $W_{[1,\alpha][1,\beta]}\sim 0$ for all in the ground state subspace which results in $F_{\text{gs}}(t\rightarrow\infty)\sim 0$. Therefore, the OTOC is dominated by the correction terms that are contributed by the excitations in the spectrum $F_{\text{ex}}(t\rightarrow\infty)$. This result is an important insight that originates from the dynamical decomposition method and cannot be observed only via real-time dynamics simulations, e.g. in Ref. Heyl et al. (2018). In addition, the operator ansatz $|W_{[1,\alpha][\theta,\beta]}|\ll 1$ guarantees a bounded correction term $F_{\text{ex}}(t\rightarrow\infty)\ll 1$. As a result, (i) the OTOC is able to capture the degeneracy in the ground state (Eq. (Detection of quantum phases via out-of-time-order correlators)) and, (ii) the correction of the excited states always remains bounded; all of which explains why the OTOC differentiates an ordered phase from a disordered one, e.g. in ground state Heyl et al. (2018) or excited-state Wang and Pérez-Bernal (2018) phase transitions. A mixed initial state (e.g. finite or infinite temperature) violates the initial state assumption, hence suggesting a smoothed phase boundary by washing away the sharp signature at the transition point Sun et al. (2018). Hence the dynamical decomposition method reveals the key ingredients of the emergent relation between information scrambling and symmetry-breaking phase transitions, rendering this unexpected numerical observation Heyl et al. (2018) a fundamental connection.

Advanced numerical methods (Lanczos, tensor networks) can be employed to determine only the lowest-lying states to give the leading order term in OTOC, Eq. (Detection of quantum phases via out-of-time-order correlators). In this sense, Eq. (Detection of quantum phases via out-of-time-order correlators) provides us a low-cost alternative to simulating the real-time OTOC dynamics in the computation of the OTOC saturation value when we use the OTOC to probe criticality. Finally, we predict that the ground state contribution to the OTOC saturation cannot efficiently distinguish quasi-long range order from a disordered phase. Because, the quasi-long range order produces zero expectation value for the order parameter (per site): $W_{[1,\alpha][1,\beta]}\sim 0$, similar to a disordered phase, and hence $F_{\text{gs}}(t\rightarrow\infty)\sim 0$ follows with correction term dominating the OTOC saturation $F(t\rightarrow\infty)$. In the following we will provide verification for our method and theory on the 1D XXZ model.

Model and results. The Hamiltonian of the XXZ model reads,

[TABLE]

where $\sigma_{i}^{n}$ are spin-1/2 Pauli matrices with energy scale set to $J$ and hence time scale set to $1/J$; $J_{z}/J$ and $h$ are the z-axis spin coupling strength and the magnetic field strength, respectively. The red lines in Fig. 1 show the phase boundaries produced by an exact method (Bethe Ansatz) for an infinite-size system. Therefore, this model has three phases: two gapped Ising phases (ferromagnetic and antiferromagnetic) at large $|J_{z}/J|$ and a gapless XY phase with quasi-long range order for small $|J_{z}/J|$, i.e. the Berezinskii-Kosterlitz-Thouless transition Berezinskiǐ (1971); Kosterlitz and Thouless (1972). We choose the OTOC operators as $\sigma^{z}_{i}$ or $\sigma^{x}_{i}$ for the spin at the $i$th site, based on the order parameters of the ferromagnetic Ising phase ($\sum_{i}\sigma^{z}_{i}$), antiferromagnetic Ising phase ($\sum_{i}(-1)^{i}\sigma^{z}_{i}$), and the XY-phase ($\sum_{i}\sigma^{x}_{i}$). Fig. 1 shows the phase diagram based on the saturation values of OTOCs with $\sigma^{z}_{i}$ [computed using Eq. (Detection of quantum phases via out-of-time-order correlators) for a system of $N=14$ spins]. We numerically confirm our theory with OTOC saturation values that are either nonzero or nearly zero in the Ising and XY phases, respectively. In fact, the OTOC recovers the phase boundaries of the Bethe ansatz solution: the agreement is perfect at the ferromagnetic-XY phase boundary and approximate at the antiferromagnetic-XY boundary due to significant finite-size effects sup .

We plot two cross-sections from Fig. 1 in Fig. 2a where the lines with orange-squares ($h/J=0$) and blue-circles ($h/J=4$) are the saturation values, Eq. (Detection of quantum phases via out-of-time-order correlators) for a short-time $tJ\sim\frac{\pi}{4}10^{1}$ (long-time results in sup ). We also plot the leading order term in the saturation, $F_{\text{gs}}(t\rightarrow\infty)$ in Fig. 2a with purple-cross ($h/J=0$) and red-diamond ($h/J=4$) lines. The OTOC saturation exactly reduces to the ground state contribution with no correction $F_{\text{ex}}=0$ in the Ising-ferromagnet, meaning that the saturation value in the ordered phase is exactly predicted by the Eq. (Detection of quantum phases via out-of-time-order correlators). The reason follows as: the spins are fully polarized in the ferromagnetic ground states, and they belong to the opposite magnetization sectors of the Hamiltonian which has magnetization conservation $[H,S_{z}]=0$ ($S_{z}=\sum_{i}\sigma^{z}_{i}$). Since they are the only states of their corresponding magnetization sectors, the fluctuations in the matrix elements are exactly zero, $|W_{[1,\alpha][\theta,\beta]}|=0$. This is why the system does not scramble information at all $F(t\rightarrow\infty)=1$, even though the XXZ model is an interacting model. We emphasize that this nonscrambling is not due to integrability of XXZ model, rather it is the signature of critical order. The rotational symmetry also protects the ferromagnetic ground states from hybridizing, all of which results in no finite-size effects at the phase boundary from ferromagnet to XY-paramagnet. In the disordered-XY phase ($h/J=0$), the ground state contribution is zero $F_{\text{gs}}=0$, leaving the correction term to dominate the saturation value, however with a small magnitude as explained above. This is the reason of the mismatch between the OTOC saturation value and its leading order term, seen in the XY-phase of Fig. 2a, while we are still able to differentiate the disordered phase from the ordered phases. Finally, in the Ising-antiferromagnet the exact agreement between Eqs. (Detection of quantum phases via out-of-time-order correlators) and (Detection of quantum phases via out-of-time-order correlators) takes place only at the $J_{z}/J\rightarrow\infty$ limit. As we approach the phase boundary towards the XY-phase, the fluctuations between matrix elements gradually increase, $|W_{[1,\alpha][1,\beta]}|\rightarrow 0$ sup , result in a nonzero but small correction term to the ground-state contribution and eventually drive the phase transition. Since the finite-size effect is significant for small sizes with exact methods, we apply density-matrix renormalization group (DMRG) algorithm with matrix product states (MPS) ITe ; sup to a system with $N=60$ and compute the real-time dynamics averaged over a short-time interval of $tJ=10$ shown with blue-circles in Fig. 2b with orange-diamonds being $F_{gs}$, Eq. (Detection of quantum phases via out-of-time-order correlators). Note that the transition point significantly shifts towards the equilibrium phase transition point, $J_{z}/J=1$. We extract the system-size scaling parameters from our DMRG computations, Fig. 2c and observe that the system indeed approaches to the equilibrium transition point when $N\rightarrow\infty$, $J_{z}^{c}=aN^{\xi}+J_{z}^{\infty}$ with exponent $\xi=-0.98$ and $J_{z}^{\infty}=1.02$ with a power-law scaling.

We plot the OTOC with $\sigma^{x}_{i}$ operator for $N=13$ (blue-circles) in Fig. 2d: the OTOC saturation remains small in all three phases and thus the OTOC can hardly distinguish the XY-ordered from XY-disordered phases. When the chains with even number of spins are used ($N=14$, orange-squares) in the theory, we do not even obtain any difference between the phases. This is in agreement with our theoretical predictions discussed previously. Additionally, the fluctuations between the matrix elements of quasi-long range order parameter $\sigma^{x}_{i}$ are always maximal regardless of the phase. Hence, we observe the mismatch between the OTOC saturation and its ground state contribution (red-diamonds $N=13$ and purple crosses $N=14$). The inset of Fig. 2d shows that the OTOC saturation value and its ground state contribution both decrease with the system size for odd-numbered chains, exhibiting that the OTOC saturation cannot capture the quasi-long range order in bigger systems and thermodynamic limit. We briefly note that the detection of the order at $J_{z}/J=-1$ is robust due to the massive degeneracy in the ground state at this point of different symmetry (SU(2) symmetry).

Conclusion. Our theoretical predictions on the XXZ model can be experimented with cold atoms Duan et al. (2003). Based on the studies in the literature Heyl et al. (2018); Wang and Pérez-Bernal (2018); Sun et al. (2018) and our results in the XXZ model, our method seems to be universal in explaining the reasoning behind the relation between scrambling and the quantum criticality. In this sense, our method is an analogue of the Eigenstate Thermalization Hypothesis: It tells us the criteria of how scrambling probes criticality; though it is independent of the integrability of the system, unlike ETH. Dynamical decomposition of OTOC is a complementary tool to the real-time evolution of a state in determining the OTOC saturation value. However in addition to providing the saturation value, it also presents us the conditions for OTOC to show either order or disorder. Based on this fact, the leading order term in our theory, Eq. (Detection of quantum phases via out-of-time-order correlators), could mark the phase transition points via system-size scalings. In conclusion, given that the initial state of OTOC is a state where the phase transition is expected to happen and the off-diagonal matrix elements of the operator are sufficiently suppressed in this state (or degenerate state subspace), OTOC could be used to dynamically detect the quantum phases with long-range order and capture the symmetry-breaking quantum phase transitions.

I Energy-time relation and finite-size effects

Eq. 5 in the main text clearly exhibits how the OTOC can suffer from the finite-size effects. Imagine that the finite-size lifts the degeneracy. Based on the dynamical equation of OTOC,

[TABLE]

one would write the ground state contribution to the OTOC at zero temperature as

[TABLE]

Since Eq. (I) is the dominant contribution to the OTOC in the ordered phase, the order will eventually be invisible to the saturation value, and it will be encoded in the frequency spectrum of the OTOC. However since this is a finite-size effect, we expect to see finite saturation value for all times in the thermodynamic limit. The period of the emerging oscillation (due to degeneracy-lifting) is $\tau=\pi/(E_{[1,2]}-E_{[1,1]})$. Then starting from $t\sim\tau/2$, the order will be invisible to the saturation value. Thus, the region where the system seems to have reached its most correlated state before the finite-size effects show up, could be defined for $t\ll\tau/2$; whereas the order will be most visible to the frequency spectrum around $t\gg\tau/2$. Further, the ground state contribution will exist in the saturation value as a non-zero effect for a time $t\sim\tau/4$, where the time-averaging will reveal the order. The relation between evolution time and energy spectrum reflects the observation that longer the time evolution, better the resolution of the energy spectrum. This, in turn, helps us to estimate the time interval of the corresponding dynamics simulation of the theory Eq. 5 (in the main text), even though Eq. 5 is explicitly time-independent.

For a phase transition that involves antiferromagnetic order, due to the doubling of the unit cell size, it is not uncommon that the finite-size contributions may oscillate as we increase the system size, depending on whether the system is composed of odd (Fig. S1) or even (Fig. 2a in the main text) number of sites. For periodic boundary conditions, the systems with even number of sites usually show smaller finite-size effects. In our studies, though, we observe the opposite: chains with odd number of spins experience finite-size effects less. As far as the OTOC is concerned, stronger finite size effects are expected for systems with even number of sites, since the key to obtain the OTOC saturation value is to sum over all the quantum states in the ground state subspace. In an Ising-ordered phase, an exact two-fold degeneracy is guaranteed for a chain with odd number of spins by the Kramers degeneracy theorem, because an odd number of spin-1/2 results in a half-integer total spin. For a system with even number of sites (i.e. integer total spin), such an exact degeneracy is not expected and thus the degeneracy is lifted by finite size effects more strongly than chains with odd number of spins. This explains the dramatic difference between the results of open and periodic boundary conditions at the XY-antiferromagnet boundary. We also note that the ground states belong to the $S_{z}=0$ magnetization sector in the antiferromagnet, which is the biggest sector of the Hamiltonian and hence they would hybridize with each other. To alleviate the finite-size effects, we make use of the time-energy relation explained above and plot the OTOC for significantly smaller interval of time, $tJ\sim\frac{\pi}{4}10$ in the main text. This should be compared with the results of a long-time evolution $tJ\sim\frac{\pi}{4}10^{3}$ in Fig. S2. This comparison is a good example of how finite-size effects could show up in scrambling, restricting the order to short-times. We also note that even though our method is valid at the infinite-time limit, due to the energy-time relation employed in the computations it produces sufficiently good results for the real-time dynamics in a short-time evolution.

Fig. S3 is the phase diagram of odd-numbered chains with phase boundaries dictated by the Bethe ansatz for an infinite-size chain. Even though the finite-size effects for small fields are more severe in even-numbered chains than odd-numbered chains in small systems, we note that the transition in odd-numbered chains is ambiguous. The continuous transition between XY- and antiferromagnetic phases due to the nature of the doubly-degenerate ground states in odd-numbered chains (as explained above), prevents a straightforward system-size scaling based on odd-numbered chains. Therefore we chose to focus on even-numbered chains in the main text. We also note that as we compute results for bigger system sizes, the apparent odd-even effect disappears. The difference between Bethe ansatz results and the OTOC phase boundary for the anti-ferromagnetic to XY phase in high fields is also due to finite-size effects. This could be seen in Figs. S4, where we compare the ground state value in OTOC with the exact phase boundaries. This means that the ground state contribution similarly suffers from the finite-size effects as well. Hence the results point to the effect of finite-size on the ground state manifold rather than the incapability of OTOC to probe the phase transition as precisely as exact results. It is an interesting question how the system size scaling results of OTOC and its ground state contribution would compare with the existing methods of determining the phase boundary, e.g. Binder ratio, fidelity measures, determining the energy gap, spatial correlation functions etc.

II Operator Ansatz Demonstrated

Here we give the additional results of XXZ model on the relation between OTOCs and phase transitions. Fig. S5 shows the difference between the OTOC saturation values, Fig. S3 and the ground state contribution in these values, Fig. S4a. The mismatch between OTOC saturation value and ground state contribution to it is clear in XY-phase, due to the fact that the correction term of the excitations is dominant in the XY-phase, as explained in the text.

Fig. S6a shows the matrix elements of the long-range order bulk operator $\sigma_{z}^{i}$ used in the study only for a ground state $|V_{1,\alpha}|^{2}$. Note how the operator’s matrix elements satisfy the operator ansatz put forward in the text: in Ising-ordered phases we observe $|W_{[1,\alpha][1,\beta]}|^{2}\gg|W_{[1,\gamma][\theta,\gamma^{\prime}]}|^{2}$ while in the Ising-disordered phase (XY-phase) $W_{[1,\alpha][1,\beta]}\sim 0$ and $|W_{[1,\alpha][\theta,\beta]}|^{2}\ll 1$ is satisfied. Fig. S6b shows the participation ratio (PR) value of the ground state in terms of spin basis. PR is defined as

[TABLE]

where $\alpha$ are eigenstates and $n$ are the reference basis. PR is a measure of fluctuations of a state in a reference basis. We see that the ferromagnetic ground states ($J_{z}/J<-1$) are more localized compared to anti-ferromagnetic ground states ($J_{z}/J>1$), because of the subspaces that they belong to under a $S_{z}$ conserving Hamiltonian. As a result, anti-ferromagnetic ground states are more susceptible to both finite-size effects (mixing in energy levels) and the effect of the rest of the terms in the Hamiltonian. This is also the reason why OTOCs are better in capturing the transition from a ferromagnet to a XY-paramagnet compared to anti-ferromagnet to XY-paramagnet. Unless $J_{z}\gg J$, the XX- and YY-coupling terms cause the Neel states to slightly couple to the other states in $S_{z}=0$ subspace.

The operator ansatz on the matrix elements can also be seen in terms of the fluctuations in a ground state. The fluctuations of the ground state can be defined as $\left(\Delta\sigma_{z}^{i}\right)^{2}=\Braket{(\sigma_{z}^{i})^{2}}-\Braket{\sigma_{z}^{i}}^{2}$ where the expectation is taken over the ground state and seen in Fig. S6c for a system size of $N=15$. Hence, we state that the fluctuations are maximized in XY-phase, causing a dominant correction term in the OTOC saturation value. The fluctuations are zero in the ferromagnetic phase and they approach to zero in the anti-ferromagnetic phase as $J_{z}/J\rightarrow\infty$. Note that this is the case because an open boundary chain with odd-number of spins have two ground states with different $S_{z}$ quantum numbers in the anti-ferromagnetic phase. This can be seen more explicitly in the matrix elements of an odd-numbered chain in Fig. S6d. The operator ansatze are satisfied as expected, however note that the main contribution comes from the diagonal elements in the Ising-ordered phases, unlike in the even-numbered chains in Fig. S6a. Finally we note that the fluctuations are always maximum for the quasi-long range operator $\sigma_{x}^{i}$.

III Comparison of real-time dynamics with theory prediction in short times

Here we share a direct comparison between real-time dynamics of OTOC in short-time and the infinite-time saturation value that is predicted by Eq. 3 in the main text, to demonstrate that the analytical framework to predict the saturation value (or the time-average) is robust in finite times. Fig. S7a shows the average of the time-signal (over a time interval of $tJ=20$ ) with blue circles and the extend of oscillations with the error bars around the blue circles. Some of these real-time OTOC signals can be seen in Fig. S7b. Even though these signals are oscillatory and show transient features, our theory could predict the average of the signals with a good accuracy as seen with red squares in Fig. S7a. Therefore, our theory is not restricted to long times. The main reason of this robustness is due to the energy-time relation we employ in the computations, hence the saturation value Eq. 3 (in the main text) changes with the interval of time-averaging even though the Eq. 3 is explicitly time-independent (see supplement Sec. I).

IV OTOC with odd number of spins and periodic boundary conditions

We present the result for odd-numbered chain if periodic boundary condition is applied, Fig. S8. The low values and fluctuations in the anti-ferromagnetic region are a sign of how OTOC is sensitive to emerging domain walls in the ground state. Domain walls are expected in the anti-ferromagnetic ground state because of the frustration in an odd-numbered periodic chain.

V The remarks on the saturation value of OTOC

Eq. 2 in main text can show why quantum chaotic spin systems should eventually decay to zero when ETH is evoked up to some approximations. When a system follows ETH, there are two criteria to satisfy: (i) $V_{\gamma\gamma^{\prime}}\ll V_{\gamma\gamma}$, where $\gamma\neq\gamma^{\prime}$, and (ii) $V_{\gamma\gamma}$ is a smooth function of energy $E_{\gamma}$, $f(E_{\gamma})$ ($V_{\gamma\gamma}$ almost do not fluctuate) Rigol et al. (2008); Srednicki (1995). In this case, we end up with $F(t\rightarrow\infty)\sim\sum_{\alpha}c_{\alpha}^{*}b_{\alpha}|V_{\alpha\alpha}|^{3}$ up to some residue due to finite-size and conservation laws Huang et al. (2019). (We assume $V=W$ for simplicity of the argument.) We can state,

[TABLE]

because $\Ket{\psi(0)}$ is drawn from a Haar-distributed ensemble (detailed in the next paragraph). Then, under the assumption of $\text{Tr}\left(V\right)=0$, $\Bra{\psi(0)}V\Ket{\psi(0)}=\sum_{\alpha}c_{\alpha}^{*}b_{\alpha}=0$. Since $V_{\gamma\gamma}$ do not fluctuate significantly via ETH’s second criteria Rigol et al. (2008) and in fact the support of distribution of $V_{\gamma\gamma}$ shrinks around the microcanonical ensemble value in the thermodynamic limit if we assume the strong form of ETH Biroli et al. (2010), $F(t\rightarrow\infty)\rightarrow 0$ for chaotic spin systems. However, in order to extract the rate of decaying to zero e.g. power-law in chaotic spin systems with conserved quantities, and/or the residue due to finite-size, one needs a more rigorous analysis where the fluctuations of the diagonal elements $V_{\gamma\gamma}$ around the smooth function $f(E)$ are included in the second assumption of ETH Huang et al. (2019).

$\text{Tr}V\sim\Braket{\psi}{V}{\psi}$ holds for a pure state $\Ket{\psi}$ randomly drawn from uniform distribution induced by Haar measure Luitz and Lev (2017). These so-called Haar-distributed random states, unlike random product states, are close to maximally entangled states Hayden et al. (2006). By being close to maximally-entangled states, Haar-distributed random pure states behave as typical states on which canonical typicality Popescu et al. (2006); Goldstein et al. (2006) could emerge. Canonical typicality is defined as the following: Consider a system $U$ with Hilbert space dimension $d_{u}$, composed of a subsystem $S$ with dimension $d_{s}$ and its environment $E$ of dimension $d_{e}$ which is significantly larger than the system itself $d_{e}\gg d_{s}$. In a system with equiprobable states, the state of the system would be $\rho_{u}=\mathbb{I}/d_{u}$. Hence the subsystem would be in a canonical state, $\Omega_{s}=\text{Tr}_{e}\rho_{u}$. Now if we take a pure typical state for the system $\Ket{\psi}$ instead of $\rho_{u}$, the subsystem state can be written as $\rho_{s}=\text{Tr}_{e}\Ket{\psi}\Bra{\psi}$. Canonical typicality asserts that $\rho_{s}\sim\Omega_{s}$ for typical states $\Ket{\psi}$. This means that typical system state $\Ket{\psi}$ is locally indistinguishable from $\rho_{u}$. As stated in Ref. Popescu et al. (2006), canonical typicality emerges because of ‘massive entanglement between the subsystem and the environment’ which is a feature of typical states. These ideas are established in Ref. Goldstein et al. (2006) and more generally in Ref. Popescu et al. (2006) via invoking Levy’s lemma Luitz and Lev (2017) for systems with constraint of energy conservation. Later studies showed that energy constraint on the system is not required for canonical typicality to emerge Reimann (2007) and it is possible to formulate the principle for mechanical observables Reimann (2007); Sugiura and Shimizu (2012). Canonical typicality for a mechanical observable reads Sugiura and Shimizu (2012); Luitz and Lev (2017)

[TABLE]

where $\epsilon$ is a small parameter. Hence, the probability that $\Braket{\psi}{V}{\psi}$ deviates from $\text{Tr}V$ decreases exponentially in the system size. Therefore, for big enough many-body systems Haar-distributed states could very well imitate an equiprobable state, $\rho_{u}=\mathbb{I}$ in the calculation of observable expectation values. This theory has been used in numerical studies to compute OTOCs at infinite-temperature Luitz and Bar Lev (2017); Daǧ and Duan (2019).

VI Density-matrix renormalization group computations

In order to compute the real-time dynamics of OTOC with MPS, we first find the ground state of the system with DMRG where we limit our computation to maximum 10 sweeps and set it to initial state. Then we time-evolve the initial state approximately Zaletel et al. (2015). Note that this approach is accurate as long as the time increment $\delta$ is small enough. Fig. S9a compares three different methods in the calculation of the OTOC at $J_{z}/J=0.5$ (XY-phase). While ED stands for exact diagonalization, single-time step means using only real time step $\delta$ and double-time means utilizing two complex time steps at a real time step in order to decrease the scaling of the error Zaletel et al. (2015). Fig. S9b shows different $\delta$ values with double-time steps in the anti-ferromagnetic phase of the XXZ model. Due to our benchmarking results, we set $\delta=0.01$ for our time evolution computations. For the calculation of the ground state contribution to OTOC as the leading-order term in the Ising-ordered phases, we compute the first three states with lowest energies via DMRG.