Finite-size scaling of out-of-time-ordered correlators at late times

TL;DR

This paper investigates how out-of-time-ordered correlators behave at late times in finite quantum many-body systems, revealing that their residual values scale differently depending on energy conservation, thus offering insights into quantum chaos.

Contribution

It demonstrates that the late-time saturation value of out-of-time-ordered correlators scales polynomially with system size when energy is conserved, contrasting with exponential decay in non-conserving systems.

Findings

Residual correlator value scales as inverse polynomial with system size when energy is conserved.

Numerical simulations support analytical predictions about late-time behavior.

Provides new understanding of quantum chaos signatures in finite systems.

Abstract

Chaotic dynamics in quantum many-body systems scrambles local information so that at late times it can no longer be accessed locally. This is reflected quantitatively in the out-of-time-ordered correlator of local operators, which is expected to decay to zero with time. However, for systems of finite size, out-of-time-ordered correlators do not decay exactly to zero and in this paper we show that the residual value can provide useful insights into the chaotic dynamics. When energy is conserved, the late-time saturation value of the out-of-time-ordered correlator of generic traceless local operators scales as an inverse polynomial in the system size. This is in contrast to the inverse exponential scaling expected for chaotic dynamics without energy conservation. We provide both analytical arguments and numerical simulations to support this conclusion.

Tables

Table 1. Table 1: Finite-size scaling of generic late-time OTOC for various types of quantum dynamics.

types of dynamics	late-time OTOC	references
Haar random unitary	$e^{-\Theta \left(n\right)}$	[18, 34]
chaotic Hamiltonian dynamics	$1/\mathrm{poly}n$	this work
many-body localization	$\Theta \left(1\right)$	[13, 7, 5, 40, 10, 4]

Equations

$$\mathrm{Re}⟨A^{†}(t)B^{†}A(t)B⟩=1-⟨[A(t),B{]}^{†}[A(t),B]⟩/2$$

$$A(t)=\sum _{(k_1,k_2,\dots ,k_n)\in \{0,x,y,z{\}}^n}{a}_{k_1{k}_2\cdots k_n}{\sigma }_{k_1}{\sigma }_{k_2}\cdots {\sigma }_{k_n}.$$

$$\sum _{k_1,k_2,\dots ,k_n}|a_{k_1{k}_2\cdots k_n}{|}^2=1.$$

$$⟨[A(t),B{]}^{†}[A(t),B]⟩=4\sum _{k_2,k_3,\dots ,k_n}|a_{y{k}_2{k}_3\cdots k_n}{|}^2+|a_{z{k}_2{k}_3\cdots k_n}{|}^2\approx 4\cdot 0.5=2.$$

$${\mathrm{OTOC}}_{\infty }(A,B,C,D):=\underset{\tau \to \infty }{\lim }\frac{1}{\tau }{\int }_0^{\tau }\mathrm{d}t⟨AB(t)CD(t)⟩,$$

$$⟨A,B,C,D{⟩}_j=(AC{)}_{jj}{B}_{jj}{D}_{jj}+A_{jj}{C}_{jj}(BD{)}_{jj}-A_{jj}{B}_{jj}{C}_{jj}{D}_{jj}.$$

$${\mathrm{OTOC}}_{\infty }(A,B,C,D)\approx \frac{1}{d}\sum _j⟨A,B,C,D{⟩}_j.$$

$${\mathrm{OTOC}}_{\infty }(A,B,A^{†},B^{†})=\frac{⟨A{A}^{†}⟩|⟨HB⟩{|}^2+⟨B{B}^{†}⟩|⟨HA⟩{|}^2}{⟨H{H}_i⟩n}+\tilde{O}(n^{-1.5}).$$

$$E_p+E_r=E_q+E_s⟹((p=q)\mathrm{and}(r=s))\mathrm{or}((p=s)\mathrm{and}(r=q)).$$

$$⟨AB(t)CD(t)⟩=\frac{1}{d}\sum _{p,q,r,s}{A}_{pq}{B}_{qr}{C}_{rs}{D}_{sp}{e}^{i(E_q-E_r+E_s-E_p)t}.$$

$${\mathrm{OTOC}}_{\infty }(A,B,C,D)=\frac{1}{d}\sum _{p,q,r,s}{A}_{pq}{B}_{qr}{C}_{rs}{D}_{sp}{\delta }_{E_p+E_r,E_q+E_s},$$

$${\mathrm{OTOC}}_{\infty }(A,B,C,D)=\frac{1}{d}\sum _{j,k}{A}_{jj}{B}_{jk}{C}_{kk}{D}_{kj}+\frac{1}{d}\sum _{j,k}{A}_{jk}{B}_{kk}{C}_{kj}{D}_{jj}-\frac{1}{d}\sum _j{A}_{jj}{B}_{jj}{C}_{jj}{D}_{jj}.$$

$${\mathrm{OTOC}}_{\infty }(H_1,H_i,H_1,H_i)=2⟨H_i^2{⟩}^2/n+O(n^{-2}).$$

$$\begin{gathered} \frac{1}{d}\sum _{j,k=1}^d(H_1{)}_{jj}(H_i{)}_{jk}(H_1{)}_{kk}(H_i{)}_{kj}=\frac{1}{d{n}^2}\sum _{j,k=1}^d{E}_j⟨j|H_i|k⟩E_k⟨k|H_i|j⟩ \\ =\frac{1}{d{n}^2}\mathrm{tr}\left(\sum _{j=1}^d|j⟩E_j⟨j|H_i\sum _{k=1}^d|k⟩E_k⟨k|H_i\right)=\frac{\mathrm{tr}(H{H}_{i}H{H}_i)}{d{n}^2}=\frac{1}{n^2}\sum _{j,k=1}^n⟨H_j{H}_i{H}_k{H}_i⟩. \end{gathered}$$

$$\frac{1}{n^2}\sum _{j=1}^n⟨H_j{H}_i{H}_j{H}_i⟩+O(n^{-2})=⟨H_i^2{⟩}^2/n+O(n^{-2})+O(n^{-2})=⟨H_i^2{⟩}^2/n+O(n^{-2}).$$

$$\frac{1}{d{n}^4}\sum _j{E}_j^4=\Theta (n^{-2}),$$

$$\frac{1}{d}\sum _j{E}_j^m=⟨H^m⟩=\Theta (n^{m/2}),\forall \mathrm{even}\mathrm{positive}\mathrm{integer}m,$$

$$\left|\frac{1}{d}\sum _j{E}_j^3\right|=|⟨H^3⟩|=O(n).$$

$$|\{j:|E_j|\ge n^{0.51}\}|/d\le O(n^{-0.01m}),\forall m>0.$$

$$\frac{1}{d}\sum _{j:|E_j|\ge n^{0.51}}{E}_j^2\le O(n^{2-0.01m}),\forall m>0.$$

$$|X_{jj}-f_X(E_j/n)|\le 1/\mathrm{poly}n$$

$$f_X(0)=\frac{1}{d}\mathrm{tr}X,$$

$$f_A^{\prime }(0)=\mathrm{tr}(HA)/\mathrm{tr}(H{H}_i),$$

$$\frac{1}{d}\sum _j|A_{jj}{|}^2=\frac{|\mathrm{tr}(HA){|}^2}{dn\mathrm{tr}(H{H}_i)}+O(n^{-2}),$$

$$\frac{1}{d}\sum _j|A_{jj}{|}^4=O(n^{-2}).$$

$$\frac{1}{d}\mathrm{tr}X=\frac{1}{d}\sum _j{X}_{jj}\approx \frac{1}{d}\sum _{j:|E_j|<n^{0.51}}{X}_{jj}\approx \frac{1}{d}\sum _{j:|E_j|<n^{0.51}}{f}_X(0)\approx \frac{1}{d}\sum _j{f}_X(0)=f_X(0),$$

$$\begin{gathered} \frac{1}{d}\mathrm{tr}(HA)=\frac{1}{d}\sum _j{E}_j{A}_{jj}\approx \frac{1}{d}\sum _{j:|E_j|<n^{0.51}}{E}_j{A}_{jj}\approx \frac{1}{d}\sum _{j:|E_j|<n^{0.51}}\frac{E_j^2}{n}{f}_A^{\prime }(0)\approx \frac{1}{d}\sum _j\frac{E_j^2}{n}{f}_A^{\prime }(0) \\ =\mathrm{tr}(H{H}_i)f_A^{\prime }(0)/d, \end{gathered}$$

$$f_A(E_j/n)=f_A(0)+f_A^{\prime }(0)E_j/n+0.5f_A^{\prime \prime }(0)E_j^2/n^2+O(|E_j{|}^3/n^3)$$

$$\frac{O(1)}{d}\sum _{j:|E_j|<n^{0.51}}\frac{|E_j{|}^3}{n^2}\le O(n^{-0.47}).$$

$$\frac{1}{d}\sum _j|A_{jj}{|}^2\approx \frac{1}{d}\sum _{j:|E_j|<n^{0.51}}|A_{jj}{|}^2\approx \frac{|f_A^{\prime }(0){|}^2}{d}\sum _{j:|E_j|<n^{0.51}}\frac{E_j^2}{n^2}\approx \frac{|f_A^{\prime }(0){|}^2}{d}\sum _j\frac{E_j^2}{n^2}=\frac{|\mathrm{tr}(HA){|}^2}{dn\mathrm{tr}(H{H}_i)},$$

$$\begin{gathered} \frac{O(1)}{d}\left|\sum _{j:|E_j|<n^{0.51}}\frac{E_j^3}{n^3}\right|+\frac{O(1)}{d}\sum _{j:|E_j|<n^{0.51}}\frac{E_j^4}{n^4}+1/\mathrm{poly}n\approx \frac{O(1)}{d}\left|\sum _j\frac{E_j^3}{n^3}\right|+\frac{O(1)}{d}\sum _j\frac{E_j^4}{n^4} \\ =O(n^{-2})+O(n^{-2})=O(n^{-2}), \end{gathered}$$

$$\frac{1}{d}\sum _j|A_{jj}{|}^4\approx \frac{1}{d}\sum _{j:|E_j|<n^{0.51}}|A_{jj}{|}^4\approx \frac{O(1)}{d}\sum _{j:|E_j|<n^{0.51}}\frac{E_j^4}{n^4}\approx \frac{O(1)}{d}\sum _j\frac{E_j^4}{n^4}=O(n^{-2}).$$

Full text

1 Introduction

Nonintegrable quantum many-body systems are expected to exhibit chaotic dynamics, which not only leads to thermalization but also scrambles local information into a nonlocal form. In the Heisenberg picture, the support of $A(t):=e^{iHt}Ae^{-iHt}$ for a local operator $A$ should grow with time under the chaotic dynamics. This growth is reflected in the noncommutativity of $A(t)$ and another local operator $B$ at a different site, which leads to the decay of the out-of-time-ordered correlator (OTOC) $\operatorname{Re}\langle A^{\dagger}(t)B^{\dagger}A(t)B\rangle$ [21, 35, 32, 36, 31, 11, 23, 33, 39, 24, 34, 8, 41, 19, 20]. Assume for simplicity that $A$ and $B$ are unitary. Then,

[TABLE]

so that when the commutator $[\cdots]$ grows, OTOC decays. The chaotic nature of the dynamics is reflected in the fast decay of OTOC away from $1$ in a relatively short time period and the approaching of OTOC to [math] at late times.

Why does chaotic dynamics lead to such decaying behavior of OTOC? While it is not possible to solve exactly the dynamics of nonintegrable systems in general, we might be able to extract some universal features, at least in certain limits. In a large class of chaotic systems without spatial locality (e.g., large-$N$ theories), OTOC at early time $t$ is given by $1-\epsilon e^{\lambda_{L}t}$, where $\epsilon$ is a small prefactor and $\lambda_{L}$ is a constant. Such an exponential deviation from the initial value is reminiscent of the so-called sensitive dependence on initial conditions (i.e., nearby points in phase space separate from each other over time at an exponential rate) in classical chaos. Thus, $\lambda_{L}$ may be interpreted as the Lyapunov exponent for quantum systems [19]. In chaotic systems with spatial locality, OTOC of two local operators starts to decay only after a delay that is proportional to the distance between the operators [32, 11, 33, 8, 41]. This is a consequence of the Lieb-Robinson bound [22, 25, 9].

In this paper, we study the behavior of OTOC at late times. For simplicity, consider a system of $n$ qubits at infinite temperature so that $\langle\cdots\rangle=\operatorname{tr}(\cdots)/2^{n}$. In the limit $t\to\infty$, a naive understanding of why OTOC approaches [math] is as follows. We expand the time-evolved operator in the $n$-qubit Pauli basis $\{\sigma_{0}=I,\sigma_{x},\sigma_{y},\sigma_{z}\}^{\otimes n}$:

[TABLE]

The unitarity of $A(t)$ implies

[TABLE]

After undergoing chaotic evolution for a sufficiently long time, the support of $A(t)$ should be the whole system, and one might expect that the coefficients $a_{k_{1}k_{2}\cdots k_{n}}$ behave like random variables due to the chaotic nature of the dynamics. If we choose $B$ to be the Pauli operator $\sigma_{x}$ of qubit $1$, then half of the terms in the expansion (2) of $A(t)$ commute with $B$ and half of them do not. Thus,

[TABLE]

The approximation step follows from Eq. (3) and the fact that we sum over half of the random variables. Substituting Eq. (4) into Eq. (1), we see that OTOC approaches [math] at late times.

Equation (2) with random coefficients is a very simple way to approximate $A(t)$ for large $t$ in chaotic systems and it is oversimplified in some respects. For example, one major difference between this approximation and the exact evolution $A(t)=e^{iHt}Ae^{-iHt}$ is that the latter preserves the spectrum of $A$ while the former does not. How does this discrepancy affect our understanding of the late-time behavior of OTOC? Is it necessary to use more refined and sophisticated approximations in order to fully capture the essence of chaotic dynamics at late times?

We focus on the scaling of late-time OTOC with system size. In finite-size systems, OTOC may converge to a small but finite value, which goes to [math] when the system size goes to infinity. One might expect this residual value to be exponentially small in the system size because we sum over an exponential number of random variables in Eq. (4). However, using a more refined approximation we show that the finite-size scaling of generic late-time OTOC should be inverse polynomial. In fact, the power-law scaling is closely related to energy conservation during the time evolution, which is not captured by simply setting the coefficients in the expansion (2) to be random.

2 Results

In this section, we introduce basic definitions and provide a summary of results.

Throughout this paper, asymptotic notations are used extensively. Let $f,g:\mathbb{R}^{+}\to\mathbb{R}^{+}$ be two positive functions. One writes $f(x)=O(g(x))$ if and only if there exist positive numbers $M,x_{0}$ such that $f(x)\leq Mg(x)$ for all $x>x_{0}$; $f(x)=\Omega(g(x))$ if and only if there exist positive numbers $M,x_{0}$ such that $f(x)\geq Mg(x)$ for all $x>x_{0}$; $f(x)=\Theta(g(x))$ if and only if there exist positive numbers $M_{1},M_{2},x_{0}$ such that $M_{1}g(x)\leq f(x)\leq M_{2}g(x)$ for all $x>x_{0}$. To simplify the notation, we use a tilde to hide a polylogarithmic factor, e.g., $\tilde{O}(f(x)):=O(f(x)\operatorname{poly}\log f(x))$.

For concreteness, consider a chain of $n$ qubits or spin-$1/2$’s with total Hilbert space dimension $d=2^{n}$ governed by a translationally invariant Hamiltonian $H=\sum_{i=1}^{n}H_{i}$, where $H_{i}$ acts on spins $i,i+1$ (nearest-neighbor interaction). While our discussion is based on a one-dimensional spin system, our results do not rely on the dimensionality of the system or the degrees of freedom being spins. A minor modification of our method leads to similar results in other settings, e.g., fermionic systems in higher dimensions. Assume without loss of generality that $\operatorname{tr}H_{i}=0$ (traceless) and $\|H_{i}\|\leq 1$ (bounded operator norm).

Let $A,B,C,D$ be local (not necessarily unitary) operators with unit operator norm. The residual value of late-time OTOC is

[TABLE]

where $\langle X\rangle:=\frac{1}{d}\operatorname{tr}X$ denotes the expectation value of an operator at infinite temperature.

Let $\{|1\rangle,|2\rangle,\ldots,|d\rangle\}$ be a complete set of eigenstates of $H$ with corresponding energies $E_{1}\leq E_{2}\leq\cdots\leq E_{d}$ in nondescending order. Let $X_{jk}=\langle j|X|k\rangle$ be the matrix element of an operator in the energy eigenbasis. Define

[TABLE]

In strongly chaotic systems, we propose the following formula for late-time OTOC:

[TABLE]

Based on this formula, we argue for

• •

${\rm OTOC}_{\infty}(A,B,A^{\dagger},B^{\dagger})$ for traceless local operators $A,B$ vanishes in the thermodynamic limit $n\to\infty$.

• •

In finite-size systems, OTOC $\langle AB(t)A^{\dagger}B^{\dagger}(t)\rangle$ saturates to $\Theta(1/n)$ if either $A$ or $B$ (or both) has a finite overlap with the Hamiltonian $H$. We not only derive the prefactor hidden in the big-Theta notation, but also provide a (not necessarily tight) upper bound on the remainder:

[TABLE]

This is our main result. It is an example where certain properties of quantum chaotic systems can be calculated analytically. For comparison, Table 1 summarizes the finite-size scaling of late-time OTOC of generic traceless local operators for various types of quantum dynamics.

The remainder of this paper is organized as follows. In Section 3, assuming a “generic” energy spectrum we present a simple derivation of Eq. (8) for the special case where the local operators in OTOC are terms in the Hamiltonian. In Section 4, we extend the approach to the general case using the eigenstate thermalization hypothesis (ETH) [6, 37, 30]. Thus, we give a rigorous proof of Eqs. (7), (8) based on two very mild assumptions for chaotic systems: a generic spectrum and ETH. In Section 5, we propose a heuristic physical picture for our results from the perspective of interpreting chaotic dynamics with random unitaries. We first introduce a previous approach, which takes into account the unitarity of the dynamics by approximating the time evolution operator $e^{-iHt}$ with a random unitary. Unfortunately, this approximation remains too crude, for it still suggests that the residual value of late-time OTOC is exponentially small in the system size. We show that once energy conservation is also taken into account by requiring the random unitary to act within small energy windows, the finite-size scaling of late-time OTOC becomes inverse polynomial. In Section 6, we support our analytical arguments with numerical simulations of a nonintegrable spin chain. The numerical results suggest that the remainder in Eq. (8) can be improved to $O(n^{-2})$.

3 Special case

In the case where the local operators in OTOC are terms in the Hamiltonian, we give a simple rigorous proof of Eq. (8) assuming only a generic spectrum.

In strongly chaotic systems, one might expect that the energy spectrum satisfies the “generic” condition:

Assumption 1** (generic spectrum; see, e.g., Ref. [38]).**

[TABLE]

This assumption is necessary in the sense that it rules out certain integrable (e.g., free-fermion) systems.

Writing out the matrix elements,

[TABLE]

Substituting into Eq. (5), we obtain

[TABLE]

where $\delta$ is the Kronecker delta. Assumption 1 implies

[TABLE]

Given a Hamiltonian $H$, there are multiple ways to write it as a sum of local terms: $H=\sum_{i}H_{i}$. Without loss of generality, we fix this ambiguity by expanding $H$ in the Pauli basis and assigning all Pauli string operators starting at site $i$ to $H_{i}$ (see Eq. (51) for an example). This convention implies $\operatorname{tr}(H_{j}H_{k})=0$ for $j\neq k$. Hence, $\langle H_{i}^{2}\rangle=\langle HH_{i}\rangle=\langle H^{2}\rangle/n$ for any $i$ due to translational invariance. Using this convention,

Theorem 1**.**

Assumption 1 implies

[TABLE]

Proof.

We will use the observation that $(H_{i})_{jj}=E_{j}/n$ for any $i$ due to translational invariance. For the present choice of local operators in OTOC, the first term on the right-hand side of Eq. (12) reads

[TABLE]

In the last sum, there are $n^{2}$ terms, most of which are zero because $\operatorname{tr}H_{j}=\operatorname{tr}H_{k}=0$. Furthermore, the convention stated above implies $\operatorname{tr}(H_{j}H_{k})=0$ for $j\neq k$. Hence, the number of nonvanishing terms in the last sum of Eq. (14) is $n+O(1)$ ($n$ comes from the terms with $j=k$ and $O(1)$ accounts for the remainder). Equation (14) equals

[TABLE]

The second term on the right-hand side of Eq. (12) gives the same result. The last term on the right-hand side of Eq. (12) equals

[TABLE]

where we used Eq. (17) with $m=4$. This completes the proof. ∎

4 General case: Implications of eigenstate thermalization

In this section, we provide an argument for Eqs. (7), (8). The argument is rigorous assuming a generic spectrum and ETH.

Lemma 1** (moments).**

[TABLE]

Proof.

Expanding $H$ in the Pauli basis, it suffices to count the number of terms that do not vanish upon taking the trace in the expansion of $H^{m}$ or $H^{3}$. ∎

Lemma 2** (concentration of eigenvalues).**

Almost all eigenstates have zero energy density:

[TABLE]

Proof.

It follows from Eq. (17) and Markov’s inequality. ∎

This lemma allows us to upper bound the total contribution of all eigenstates away from the middle of the spectrum, e.g.,

[TABLE]

Lemmas 2 and Eqs. (17), (20) are related to the fact that $E_{j}$’s approach a normal distribution in the thermodynamic limit $n\to\infty$ [14, 3]. Indeed, $|E_{j}|=\Theta(\sqrt{n})$ for almost all $j$.

It suffices to assume ETH for eigenstates in the middle of the spectrum.

Assumption 2** (eigenstate thermalization hypothesis in the middle of the spectrum).**

Let $\delta$ be an arbitrarily small positive constant. For any local operator $X$ with $\|X\|\leq 1$, there is a function $f_{X}:[-\delta,\delta]\to[-1,1]$ such that

[TABLE]

for all $j$ with $|E_{j}|\leq\delta n$, where $\operatorname{poly}n$ denotes a polynomial of sufficiently high degree in $n$. We assume that $f_{X}$ is smooth in the sense of having a Taylor expansion to some low order.

It was proposed analytically [38] and supported by numerical simulations [17] that the right-hand side of Eq. (21) can be improved to $e^{-\Omega(n)}$. For our purposes, however, a (much weaker) inverse polynomial upper bound suffices.

Lemma 3**.**

For any local operator $X$ and traceless local operator $A$, Assumption 2 implies

[TABLE]

For a generic traceless local operator $A$, the right-hand side of Eq. (23) (the normalized overlap between $A$ and the Hamiltonian) is finite and the first term on the right-hand side of Eq. (24) is $\Theta(1/n)$.

Proof of Eq. (22).

[TABLE]

where we used Lemma 2 in the second and fourth steps. The third step follows from the continuity of $f_{X}(x)$ at $x=0$. Taking the limit $n\to\infty$, all errors in Eq. (26) vanish and thus we obtain Eq. (22). In particular, $f_{A}(0)=0$ for any traceless local operator $A$. ∎

Proof of Eq. (23).

[TABLE]

where we used Lemma 2 and Eq. (20) in the second and fourth steps, respectively. In the third step, we used Eq. (21) and the Taylor expansion

[TABLE]

so that the approximation error in this step is upper bounded by

[TABLE]

Taking the limit $n\to\infty$, all errors in Eq. (27) vanish and thus we obtain Eq. (23). ∎

Proof of Eq. (24).

[TABLE]

where we used Lemma 2 and Eqs. (20), (23) in the first, third, and last steps, respectively. In the second step, we used Eqs. (21), (28) with the approximation error upper bounded by

[TABLE]

where we used Eqs. (17), (18). ∎

Proof of Eq. (25).

[TABLE]

∎

Let $J\subseteq\mathbb{R}$ be an energy interval. Define

[TABLE]

as the projector onto $J$.

Lemma 4** ([1]).**

Let $\epsilon<\epsilon^{\prime}$. For any local operator $X$,

[TABLE]

This lemma states that local operators cannot (up to an exponentially small error) connect projectors that are far away from each other in the spectrum.

Justification of Eq. (7).

Let $c$ be a sufficiently large constant. Consider the first term on the right-hand side of Eq. (12):

[TABLE]

where we used Lemma 4 in the first and third steps: Due to the presence of off-diagonal matrix elements $B_{jk},D_{kj}$, the total contribution of all terms with $|E_{j}-E_{k}|\geq c\ln n$ is upper bounded by $1/\operatorname{poly}n$. In the second step of Eq. (35), we replace $C_{kk}$ by $C_{jj}$ using ETH (Assumption 2), which states that eigenstates with similar energies have similar local expectation values. A detailed and rigorous error analysis for Eq. (35) is given in Propositions 1, 2 below.

Equation (35) shows that the first term on the right-hand side of Eq. (12) corresponds to the second term on the right-hand side of Eq. (6). Similarly, the second term on the right-hand side of Eq. (12) corresponds to the first term on the right-hand side of Eq. (6). Obviously, the last terms on the right-hand sides of Eqs. (6), (12) are the same. Thus, we obtain Eq. (7). ∎

Proposition 1**.**

The approximation errors in the first and third steps of Eq. (35) are $1/\operatorname{poly}n$, where $\operatorname{poly}n$ denotes a polynomial of sufficiently high degree in $n$.

Proof.

Let

[TABLE]

Since $\tilde{C}$ is the diagonal part of $C$ (in the energy eigenbasis), it is easy to see $\|\tilde{C}\|\leq\|C\|$. The approximation error in the first step of Eq. (35) is

[TABLE]

where we used Lemma 4. The approximation error in the third step of Eq. (35) can be upper bounded similarly. ∎

Proposition 2**.**

The approximation error in the second step of Eq. (35) is $\tilde{O}(n^{-1.5})$.

Proof.

Let $n$ be sufficiently large such that $n^{0.51}+c\ln n<\delta n$, and define

[TABLE]

For $j,k$ such that $|E_{j}|<n^{0.51}$ and $|E_{j}-E_{k}|<c\ln n$, Assumption 2 implies

[TABLE]

Hence, $\|\tilde{C}^{(j)}\|=\tilde{O}(1/n)$ for any $j$ such that $|E_{j}|<n^{0.51}$. The approximation error in the second step of Eq. (35) is

[TABLE]

where we used Eq. (24) in the last step. ∎

Justification of Eq. (8).

Specializing to $\langle AB(t)A^{\dagger}B^{\dagger}(t)\rangle$, the derivation above yields

[TABLE]

Consider the first term on the right-hand side:

[TABLE]

where we used Lemma 2 in the first and third steps; the continuity of $f_{AA^{\dagger}}(x)$ at $x=0$ in the second step; Eqs. (22), (24) in the last step. A rigorous error analysis for Eq. (42) is given in Proposition 3 below.

The second term on the right-hand side of Eq. (41) can be estimated similarly. The third term on the right-hand side of Eq. (41) is

[TABLE]

where we used Eq. (25). Thus, Eq. (8) is proved based on Assumptions 1, 2. ∎

Proposition 3**.**

The error in Eq. (42) is $O(n^{-2})$.

Proof.

The approximation error in the last step of Eq. (42) is $O(n^{-2})$ as given by Eq. (24). Using the Taylor expansion of $f_{AA^{\dagger}}(x)$ at $x=0$, we estimate the approximation error in the second step of Eq. (42):

[TABLE]

where we used the Taylor expansion of $f_{B}(x)$ at $x=0$ and the inequality of arithmetic and geometric means in the first step; Lemma 2 in the second step; Eqs. (17), (18), (25) in the third step. ∎

5 Chaotic dynamics as random unitary

In this section, we rederive Eq. (7) using techniques from the theory of random unitaries. The derivation is not rigorous, but provides a heuristic picture showing the extent to which chaotic dynamics can be approximated by a random unitary.

To improve the approximation described by Eq. (2), we first take into account the unitarity of the dynamics. In strongly chaotic systems, it is tempting to expect

Assumption 3**.**

The time evolution operator $e^{-iHt}$ for large $t$ behaves like a random unitary.

Based on this assumption, late-time OTOC can be estimated from

[TABLE]

where $U$ is taken from the unitary group $\mathcal{U}(d)$ with respect to the Haar measure.

Lemma 5** ([18, 34]).**

[TABLE]

where $\langle XY\rangle_{c}:=\langle XY\rangle-\langle X\rangle\langle Y\rangle$ is the connected correlator and

[TABLE]

Note that the right-hand side of Eq. (6) resembles that of Eq. (47) in the sense of replacing every $\langle\cdots\rangle$ (expectation value at infinite temperature) by $\langle j|\cdots|j\rangle$ (expectation value in an eigenstate).

Corollary 1** ([18, 34]).**

Assumption 3 and Lemma 5 imply

[TABLE]

Therefore,

• •

${\rm OTOC}_{\infty}(A,B,A^{\dagger},B^{\dagger})$* for traceless operators $A,B$ vanishes in the thermodynamic limit $n\to\infty$.*

• •

In finite-size systems, the saturation value of OTOC $\langle AB(t)A^{\dagger}B^{\dagger}(t)\rangle$ is exponentially small in the system size (because $d=2^{n}$).

The approximation stated in Assumption 3 is still too crude. We propose a refined version of Assumption 3 by incorporating energy conservation and argue (nonrigorously) that Eq. (7) follows from this refinement.

We observe that the time evolution conserves energy and that local operators can only additively change the energy of a state by $O(1)$ (Lemma 4). Thus, the action of OTOC $AB(t)CD(t)$ is approximately restricted to each microcanonical ensemble. This observation motivates a refinement of Assumption 3 in strongly chaotic systems:

Assumption 4**.**

The time evolution operator $e^{-iHt}$ for large $t$ behaves like a random unitary in each microcanonical ensemble.

Conceptually, this assumption is related to the so-called random diagonal unitaries [27, 26].

Based on Assumption 4, we argue for Eq. (7). Since the bandwidth of $H$ is $\Theta(n)$, we decompose the energy spectrum into a disjoint union of $\Theta(n/\Delta)$ microcanonical ensembles with bandwidth $\Delta$. Let $J_{k}:=[k\Delta,(k+1)\Delta)$ and define $[A,B,C,D]_{k}$ as the right-hand side of Eq. (47) with every $\langle\cdots\rangle$ replaced by the expectation value $\operatorname{tr}(P_{J_{k}}\cdots)/\operatorname{tr}P_{J_{k}}$ in the microcanonical ensemble. We expect

[TABLE]

The first step is a consequence of Lemma 5 and Assumption 4. Indeed, it is just Eq. (48) restricted to the microcanonical ensemble $P_{J_{k}}$. The last step of Eq. (49) used ETH. Equation (7) follows immediately from Eq. (49).

An important subtlety here, which does not appear in the derivation of Eq. (48), requires further explanation. For an eigenstate $|j\rangle$ in a microcanonical ensemble $P_{J_{k}}$, the state $AB(t)CD(t)|j\rangle$ may not be completely in the microcanonical ensemble. As long as $A,B,C,D$ are local operators, Lemma 4 implies

[TABLE]

i.e., the “leakage” out of the microcanonical ensemble $P_{J_{k}}$ is exponentially small. This is why Eq. (7) requires the locality of $A,B,C,D$, although Corollary 1 does not.

6 Numerics

In this section, we support Eq. (8) with numerical simulations. Consider the spin-$1/2$ chain

[TABLE]

with periodic boundary conditions ($\sigma_{n+1}^{z}:=\sigma_{1}^{z}$), where $\sigma_{i}^{x},\sigma_{i}^{y},\sigma_{i}^{z}$ are the Pauli matrices at site $i$. For $g=0$, this model is nonintegrable in the sense of Wigner-Dyson level statistics [2, 16]. Reference [32] calculated OTOC, focusing on the butterfly effect rather than the late-time behavior. Note that for $g=0$, most energy levels are twofold degenerate so that Assumption 1 does not hold.

We fix $g=0.1$. Intuitively, the model is nonintegrable for any value of $g$. We have numerically confirmed the validity of Assumption 1 for $n=5,6,\ldots,12$. Presumably, Assumption 1 holds for any integer $n\geq 5$. Let

[TABLE]

Note that the values of $F_{n}^{x},F_{n}^{z}$ are independent of $i$. We compute $F_{n}^{x},F_{n}^{z}$ using exact diagonalization. The results are shown in the left panel of Fig. 1.

The leading terms in the finite-size scaling of $F_{n}^{x},F_{n}^{z}$ are calculated analytically from Eq. (8):

[TABLE]

We expect that the noticeable differences between $G_{n}^{x},G_{n}^{z}$ and the power-law fits to $F_{n}^{x},F_{n}^{z}$ are due to finite-size effects. To justify this claim, we perform a scaling analysis of the errors $|F_{n}^{x}-G_{n}^{x}|$, $|F_{n}^{z}-G_{n}^{z}|$ in the right panel of Fig. 1. The numerics suggest that the errors should vanish as $\Theta(n^{-2})$ in the thermodynamic limit $n\to\infty$.

7 Conclusion

We propose that in order to better approximate the late-time behavior of chaotic dynamics generated by a time-independent Hamiltonian, one needs to take into account energy conservation. In particular, we show that approximation schemes with and without energy conservation make different predictions about OTOC at late times: without energy conservation, late-time OTOC scales inverse exponentially with system size; with energy conservation, the scaling is inverse polynomial. The latter prediction has been rigorously confirmed based on two very mild assumptions and is consistent with numerical simulations of a nonintegrable spin chain.

An immediate open question is how good the energy-preserving approximation scheme proposed in this paper is in predicting the late-time behavior of higher-order time-ordered or out-of-time-ordered correlators. A more general problem for future study is how to approximate the time evolution process and capture other universal features of chaotic dynamics. See Refs. [15, 28, 29, 12] for recent progress in this direction.

Acknowledgments

We would like to thank Xie Chen, Yingfei Gu, Nick Hunter-Jones, Alexei Y. Kitaev, and Douglas Stanford for helpful discussions. We are especially grateful to X.C., who wrote a substantial portion of this paper. We acknowledge funding provided by the Institute for Quantum Information and Matter, an NSF Physics Frontiers Center (NSF Grant PHY-1733907). Additional funding support was provided by NSF DMR-1654340 (Y.H., Y.-L.Z.).