Information Scrambling and Loschmidt Echo

TL;DR

This paper establishes an analytical and numerical connection between the out-of-time-order correlator, a measure of quantum information scrambling, and the Loschmidt echo, a classical diagnostic of chaos, enhancing experimental accessibility.

Contribution

It reveals a direct relationship between the out-of-time-order correlator and the Loschmidt echo, linking quantum information scrambling to a well-known dynamical chaos measure.

Findings

Analytical demonstration of the correlator-echo relationship

Numerical verification confirming the theoretical link

Provides a practical approach for experimental studies of quantum chaos

Abstract

We demonstrate analytically and verify numerically that the out-of-time order correlator is given by the thermal average of Loschmidt echo signals. This provides a direct link between the out-of-time-order correlator -- a recently suggested measure of information scrambling in quantum chaotic systems -- and the Loschmidt echo, a well-appreciated familiar diagnostic that captures the dynamical aspect of chaotic behavior in the time domain, and is accessible to experimental studies.

Equations

$$M(t)\equiv |⟨\psi |e^{i{H}_{0}t}{e}^{-i(H_0+V)t}|\psi ⟩{|}^2.$$

$$F_{\beta }(t)\equiv ⟨A^{†}(t)B^{†}A(t)B{⟩}_{\beta },$$

$$\overline{F_{\beta }(t)}\equiv \int dAdB{F}_{\beta }(t),$$

$$\int dA{A}^{†}OA=\frac{1}{d_A}{\mathbb{I}}_A\otimes {\mathrm{tr}}_A(O),$$

$$\int dA{F}_{\beta =0}(t)\equiv \frac{1}{d}\int dA\mathrm{tr}(A^{†}(t)B^{†}A(t)B)$$

$$=$$

$$H={\mathbb{I}}_A\otimes H_B+H_A\otimes {\mathbb{I}}_B+H^{\prime },$$

$$H^{\prime }\equiv \delta \sum _k{V}_A^k\otimes V_B^k.$$

$$B(-t)$$

$$\propto \overline{e^{-i(H_B+V_{\alpha })t}B{e}^{i(H_B+V_{\alpha })t}}.$$

$$V_{\alpha }=\sum _k\pm \delta V_B^k.$$

$$\overline{F_{\beta =0}(t)}\approx \overline{|⟨e^{i(H_B+V_{\alpha })t}{e}^{-i(H_B+V_{{\alpha }^{\prime }})t}{⟩}_{\beta =0}{|}^2}.$$

$$\overline{F_{\beta =0}(t)}\approx |⟨e^{i(H_B+V_1)t}{e}^{-i(H_B+V_2)t}{⟩}_{\beta =0}{|}^2.$$

$$\sum _{i=1,2}\left(\frac{1}{2m_i}{\hat{p}}_i^2-\frac{m_i{\omega }_i^2}{2}{\hat{x}}_i^2\right)+\delta {\hat{x}}_1{\hat{x}}_2.$$

$$M(t)=\frac{1}{2}+\frac{1}{2}|⟨{\psi }_2|e^{i(H_2+\delta c_1{\hat{x}}_2)t}{e}^{-i(H_2-\delta c_1{\hat{x}}_2)t}|{\psi }_2⟩{|}^2.$$

$$M_1(t)=1-{\delta }^2\frac{4c_1^2{c}_2^2}{(i{\omega }_2{)}^2}{\sin }^2(i{\omega }_{2}t),$$

$$\int dA{A}^{†}OA=\frac{1}{d}\mathrm{tr}(O)\mathbb{I}.$$

$$\int d\mu (A)A^{†}OA=U^{†}\left(\int d\mu (A)A^{†}OA\right)U.$$

$$\mathrm{tr}\int d\mu (A)A^{†}OA=\int d\mu (A)\mathrm{tr}(O).$$

$$\int dA{A}^{†}\otimes {\mathbb{I}}_B{O}_{AB}A\otimes {\mathbb{I}}_B$$

$$=\int dA{A}^{†}\otimes {\mathbb{I}}_B\sum _i{O}_i^A\otimes O_i^{B}A\otimes {\mathbb{I}}_B$$

$$=\sum _i\int dA{A}^{†}{O}_i^{A}A\otimes O_i^B$$

$$=\frac{1}{d_A}\sum _i\mathrm{tr}(O_i^A){\mathbb{I}}_A\otimes O_i^B=\frac{1}{d_A}{\mathbb{I}}_A\otimes {\mathrm{tr}}_A{O}_{AB}.$$

$$H={\mathbb{I}}_A\otimes H_B+H_A\otimes {\mathbb{I}}_B+H^{\prime },$$

$$H^{\prime }\equiv \delta \sum _{k=1}^{d_A^2}{V}_A^k\otimes V_B^k.$$

$$\mathrm{tr}(V_A^i{V}_A^j)=d_A{\delta }_{i,j}.$$

$$B(-t)={\mathrm{tr}}_A\left(e^{-iHt}{\mathbb{I}}_A\otimes B{e}^{iHt}\right).$$

$$\delta \mathcal{F}(t)=\delta \sum _k{l}_k(t)V_B^k,$$

$$\overline{l_i(t)l_j(t-\tau )}$$

$$=$$

$$\approx$$

$$B(t)=d_A\times \overline{e^{-i(H_B+\mathcal{F})t}B{e}^{i(H_B+\mathcal{F})t}},$$

$$B(-t)\approx d_A\times \frac{1}{N}\sum _{\alpha =1}^N{e}^{-i(H_B+\delta V_{\alpha })t}B{e}^{i(H_B+\delta V_{\alpha }))t},$$

$$\int dA{F}_{\beta =0}(t)\equiv \frac{1}{d}\int dA\mathrm{tr}(A^{†}(t)B^{†}A(t)B)$$

$$=$$

Full text

Introduction.— The study of the quantum version of classically chaotic systems gave rise to the field of quantum chaos, since it was realized that quantum chaotic systems share certain common characteristics Haake (2006); Berry (1989). In particular, in spite of the absence of the telltale exponential sensitivity to initial conditions in unitary quantum evolution, one can use as a quantum diagnostic of chaos sensitivity of the evolution to small perturbations of the Hamiltonian Peres (2002) or its entropy production in the presence of coupling to the environment Zurek and Paz (1994); Zurek (1998). These and related manifestations of quantum chaos have been by now widely studied Levstein et al. (1998); Jalabert and Pastawski (2001); Cucchietti et al. (2003); Kohler and Recher (2012); Sánchez et al. (2016); Chenu et al. (2018, 2019); Sánchez et al. (2020), using diagnostics such as Loschmidt echo (LE) Goussev et al. (2012); Gorin et al. (2006),

[TABLE]

This quantity incorporates the simple idea Peres (2002) that small perturbations of the Hamiltonian may trigger dramatic changes of the dynamics, inducing the butterfly effect.

More recently, the out-of-time-order correlator (OTOC) Larkin and Ovchinnikov (1969); Kitaev , another diagnostic for quantum chaos, has been proposed and received considerable attention across many different fields in physics, including quantum information and high-energy and condensed matter physics Shenker and Stanford (2014a, b); Roberts et al. (2015); Roberts and Stanford (2015); Blake (2016); Roberts and Swingle (2016); Maldacena et al. (2016); Fu and Sachdev (2016); Hosur et al. (2016); Tsuji et al. (2018); Swingle and Chowdhury (2017); Campisi and Goold (2017); Chen et al. (2017); Rozenbaum et al. (2017); Dóra and Moessner (2017); Hashimoto et al. (2017); Lin and Motrunich (2018a, b); von Keyserlingk et al. (2018); Nahum et al. (2018); Khemani et al. (2018); Gharibyan et al. (2018); Xu and Swingle (2019); Tuziemski (2019). The OTOC is defined as a four-point correlator with unusual time ordering:

[TABLE]

where $A$ and $B$ are typically chosen as local operators; $A(t)=e^{iHt}Ae^{-iHt}$ is the Heisenberg operator evolving under total Hamiltonian $H$ and the average $\langle\rangle$ is taken over a thermal state at the inverse temperature $\beta$. In chaotic systems, the OTOC exhibits fast decay and converges to a persistent small value Kitaev . It was argued that under certain natural assumptions, the exponential decay rate is bounded by $\lambda\leq 2\pi/\beta$ Maldacena et al. (2016); Tsuji et al. (2018). Another benefit of OTOC is that it is designed to probe the spreading of local information over the entire system. Moreover, for systems with spatial structures, information measured by the OTOC propagates ballistically with a finite velocity known as the butterfly velocity Roberts and Swingle (2016); Nahum et al. (2018); Khemani et al. (2018).

The OTOC is typically understood as an intrinsic echo-type quantity. For instance, when $A$ and $B$ are chosen as unitary operators, Eq. (2) can be directly measured by echo experiments Swingle et al. (2016); Gärttner et al. (2017); Li et al. (2017); Landsman et al. (2019); Sánchez et al. (2020). There are efforts to build more direct links between these two quantities, e.g., using variants of the OTOC and LE or particular choices of operators for the OTOC evaluation Kurchan (2018); Schmitt et al. (2019); Hamazaki et al. (2018). However, the precise quantitative equivalence between the OTOC and the Loschmidt echo is still missing. Establishing such a relation would be beneficial for both areas and shed new light on the whole field of quantum chaos.

In this work, we accomplish the task of connecting the OTOC to the Loschmidt echo. We shall focus on the temporal decay of the OTOC without extra complications caused by spatial propagation. We demonstrate that the OTOC equals the thermal average of the Loschmidt echo. The coupling between the target local systems, i.e., the supports of the local operators, and the rest of the total system plays the role of a perturbation. To further support our theory, we present two model studies involving coupled inverted harmonic oscillators and a random matrix model.

*Bridging the out-of-time-order correlator and Loschmidt echo.—*For a chaotic Hamiltonian, the universal decay of the OTOC is insensitive to the form of operators $A$ and $B$ in Eq. (2), as long as they are generic, i.e., not reflecting the particular symmetries possessed by the Hamiltonian. Any generic choice of local operators, e.g., random operators, is representative for the universal decay of the OTOC. This allows us to look at the typical behavior of the OTOC by averaging all unitary operators on subsystem $S_{A}$ and $S_{B}$:

[TABLE]

where the integration is performed with respect to the Haar measure for unitary operators. Similar ideas have been considered in the literature Fan et al. (2017); Cotler et al. (2017a); de Mello Koch et al. (2019); Ma (2019). Here, we will assume that $A$ and $B$ are supported on distinct local subsystems. For global operators the OTOC has been shown to be closely related to the spectral form factor of the Hamiltonian Cotler et al. (2017a, b); de Mello Koch et al. (2019). As will be seen in the following, taking into account the local structure of the system is crucial to reveal the correct behavior of the OTOC.

For simplicity, we focus on the OTOC at an infinite temperature ($\beta=0$). It is straightforward to generalize to a finite temperature by distributing the thermal density operator over a thermal loop, e.g., using the scheme described in Refs. Maldacena et al. (2016); Cotler et al. (2017a); Tsuji et al. (2018). The finite-temperature correction will be taken into account when discussing the temperature dependence of the decay rate.

We focus on the scenario that $A$ is an operator with support on a small local subsystem $S_{A}$, while operator $B$ is chosen such that its support $S_{B}$ is the complement of $S_{A}$, as illustrated in Fig. 1. It is reasonable to expect that choosing operator $B$ in such a manner also captures the spreading of operator $A$ over the entire system, detected by its nonzero projection at later times on the support of operator $B$. Analysis of this particular scenario is also instructive and can be generalized in a similar way to cases where $B$ is a small local operator as well.

The Haar integral over subsystem operators can be evaluated with the aid of the formula

[TABLE]

where $\mathbb{I}_{A}$ is the identity operator and $\operatorname{tr}_{A}$ is the partial trace over subsystem $S_{A}$. The Haar measure is unique up to a constant multiplicative factor. Here, $d_{A}$, the dimension of $S_{A}$ is introduced as a convention to normalize the OTOC Bhattacharyya et al. (2019). The proof of the above equation is given in Appendix A.

The average over all random unitary operators on subsystem $S_{A}$ gives us

[TABLE]

The last line of the above equation involves the reduced dynamics of operator $B$, i.e., $B(-t)\equiv\operatorname{tr}_{A}(e^{-iHt}Be^{iHt})$. In order to further perform the average over $B$, we estimate $B(-t)$ in the following way. The total system Hamiltonian, in general, has the structure

[TABLE]

Here $V^{k}_{A}$’s are Hermitian and orthonormal (with respect to the Hilbert-Schmidt inner product and norm); $V^{k}_{B}$’s are Hermitian, orthogonal and have norms on the same order as $H_{B}$, such that $\delta$ quantifies the relative strength of the coupling. In realistic physical systems, the coupling parameter $\delta\ll 1$. To get the solution of $B(-t)$, we replace the effect of the coupling with an ensemble of random noises $\{V_{\alpha}\}$ on subsystem $S_{B}$, namely,

[TABLE]

Here $\alpha$ labels different realizations of the noise and the average taken is over all realizations of the noise operator $V_{\alpha}$’s.

The above claim is based on the correspondence between the symptoms of decoherence Zurek (2003); Schlosshauer (2007) (process that involves entangling correlations between the system and the environment) and symptoms of the suitable external noises. In Appendix B, we back up this claim with a more mathematically rigorous treatment of the noise operators in terms of master equations up to the second order of the coupling parameter $\delta$ Budini et al. (1999), where it is shown that the noise operators can be chosen as linear combinations of $V_{B}^{k}$’s with random coefficients $\pm 1$:

[TABLE]

With the aid of the alternative form for the reduced dynamics for $B(-t)$ in Eq. (7), averaging over operators $B$ can be further performed in the same manner as for the operators $A$. This gives the final expression for the averaged OTOC:

[TABLE]

Equation (9) is the main result of this work. Note that: (i) As has been mentioned before, the above result generalizes to finite temperatures. We present the full derivation in Appendix C, using the same scheme for regularizing the thermal state explored in Refs Maldacena et al. (2016); Cotler et al. (2017a); Tsuji et al. (2018). This regularization scheme is crucial for the discussion of the bound of the OTOC decay Romero-Bermúdez et al. (2019). (ii) The perturbations that appear in the LE emerge from the interactions nested in the total Hamiltonian. (iii) The OTOC has several decay regimes, e.g., the early growth $\sim a-\epsilon be^{\lambda t}+O(\epsilon^{2})$ before the scrambling (Ehrenfest) timescale Maldacena et al. (2016), where $a$ and $b$ are order-1 numbers and $\epsilon\ll 1$ is a small parameter; or the intermediate pure exponential decay. The only approximation involved in the derivation of Eq. (9) is the second-order approximation of the coupling parameter $\delta$. It will be shown in the following that the OTOC-LE connection holds in both the scrambling and intermediate decay regimes. Remarkably, the second order of the coupling parameter $\delta^{2}$ is identified as the small prefactor $\epsilon$ in the early exponential growth, which determines the scrambling timescale $\sim(1/\lambda)\ln(1/\epsilon)$.

Case I: Early scrambling.— To verify the OTOC-LE relation in the scrambling regime, we first study an exactly solvable model consisting of two coupled inverted harmonic oscillators (IHOs), whose Hamiltonian is

[TABLE]

This model employs IHOs, as they were first used to emulate dynamical instability characteristic of quantum chaos in an exactly solvable system Blume-Kohout and Zurek (2003). Their role here is to capture the essential ingredients of the OTOC-LE relation: The OTOC is computed for two (local) operators on the two oscillators, respectively. The parameters are tuned as $m_{1}\gg m_{2}$, $\omega_{1}\ll\omega_{2}$ and $\delta\ll 1$ such that the dynamics of the first oscillator is much slower than the second one and the coupling is weak. Hence, oscillator $1$ and $2$ mimic the subsystems $S_{A}$ and $S_{B}$, respectively..

Since the spectrum of the IHOs is not bounded from below, the thermal state is not well defined. We replace the regularized thermal state with a pure state $|\psi_{1}(x_{1})\rangle|\psi_{2}(x_{2})\rangle$. The pure state average of the OTOC has been considered in the literature Gärttner et al. (2017). In this case, the Haar averaged OTOC of the two coupled oscillators can be computed exactly. We depict the result in Fig. 2 and delegate the lengthy calculation to Appendix D, where it is shown that the pure state average carries out in the same manner as a thermal ensemble.

According to Eq. (9), the OTOC equals the LE of the second oscillator (the large subsystem $S_{B}$), which admits an exact solution as well. The perturbations emerging from the coupling are $\pm\delta c_{1}\hat{x}_{2}$ with equal probability, where $c^{2}_{1}=\langle\psi_{1}|x^{2}_{1}|\psi_{1}\rangle$. Note that the average of $x^{2}_{1}$ appears because of the normalization condition in the decomposition of the interaction in Eq. (6). See details about the derivation of the perturbation operators in Appendix D. As a simple system, the coupled IHO has only one noise operator. Thus the coarse-grained version of the OTOC-LE connection in Eq. (9b) is not reliable. We instead explicitly use the exact form in Eq. (9a), the right-hand side of which reduces to

[TABLE]

Denote $M_{1}(t)\equiv|\langle\psi|e^{i(H_{2}+\delta c_{1}\hat{x}_{2})t}e^{-i(H_{2}-\delta c_{1}\hat{x}_{2})t}|\psi\rangle|^{2}$. To extract the early exponential growth, consider the quantity $1-\overline{F(t)}\approx 1-M(t)=\frac{1}{2}\left(1-M_{1}(t)\right)$. This implies that, in the early growth regime, the averaged OTOC $\overline{F(t)}$ has the same exponential growth rate as $M_{1}(t)$, but a prefactor half of the latter. Note that the growth of the LE $M(t)$ does not saturate to one. The reason is that the OTOC-LE connection is exact up to second order of the coupling parameter $\delta$. As will be demonstrated in the following, $\delta^{2}$ plays the role of the prefactor $\epsilon$ of the OTOC early growth $\sim\epsilon e^{\lambda t}$, which is precisely reflected in the early growth of the LE. The saturation of the OTOC is induced by higher orders of $\epsilon$.

To see the significance of the coupling strength, we expand $M_{1}(t)$ to second order of $\delta$ using the Baker-Campbell-Hausdorff (BCH) formula, which gives

[TABLE]

where $c^{2}_{2}=\langle\psi_{2}|x_{2}^{2}|\psi_{2}\rangle$. This describes an exponential growth with rate $2\omega_{2}$. The second order of the coupling parameter plays the role of the prefactor in the early exponential Lyapunov growth. Derivations for the exact solution of $M_{1}(t)$ and its second-order BCH expansion are presented in Appendix D. Figure 2 depicts the exact solutions of the Haar averaged OTOC, $M_{1}(t)$ as well as the exponential growth extracted from the second-order expansion in Eq. (12). The OTOC-LE connection and the early scrambling are clearly revealed.

Case II: Intermediate decay.— The early exponential growth of the OTOC has been predicted and observed in various platforms (See discussions in Ref. Xu and Swingle (2020) and the references therein). However, in a variety of systems, especially finite-size lattice systems, the scrambling time is typically too small to be reliably visible. Instead, the decay of the OTOC is a pure exponential. To reveal the OTOC-LE relation in such an intermediate decay regime, we study a random matrix model. The model Hamiltonian takes the general form Eq. (6). The noninteracting part of the Hamiltonian are random matrices from the standard Gaussian unitary ensemble (GUE), whose matrix elements have independent real and imaginary parts as Gaussian random numbers with zero mean and unit variance. The small subsystem $S_{A}$ is chosen as a single qubit. The coupling matrix $V^{k}_{A}$ in the decomposition Eq. (6) are chosen as Pauli matrices, i.e., $\{I,\sigma_{x},\sigma_{y},\sigma_{z}\}$, while $V^{k}_{B}$ on the large subsystem $S_{B}$ are drawn from the GUE.

The random matrix model has enough complexity to make the OTOC insensitive to the choice of operators. This allows us to numerically simulate the evolution of the OTOC for two random Hermitian operators on the corresponding subsystems. The LE can be computed using the coarse-grained version of the OTOC-LE relation in Eq. (9b). The effective perturbation operators, by Eq. (8), are the sum of four random matrices from the GUE and, therefore, have matrix elements with variance $4\delta^{2}$.

It is well known that in the intermediate regime the decay of the LE depends on the strength of the perturbation; i.e., for small perturbations the decay is exponential, while for large perturbations the decay can be Gaussian. Since the relative coupling strength decreases with the system size, exponential decays are typically expected in the thermodynamics limit. However, for small size systems Gaussian decays might be observed. This explains the Gaussian decay of the OTOC in finite-size systems, e.g., the Sachdev-Ye-Kitaev (SYK) model Sachdev and Ye (1993); Kitaev ; Maldacena and Stanford (2016), for which in the large-$N$ limit the OTOC switches from an early growth to an intermediate exponential decay. As another application of our theory, we present a detailed discussion of the SYK model in Appendix E.

The numerical simulations for the random matrix model are presented in Fig. 3, where the OTOC and LE are shown to match very well; both exponential and Gaussian decays are revealed.

The bound on chaos.— In the scrambling regime, the early exponential growth rate was conjectured to be universally bounded by the temperature, $\lambda\leq 2\pi/\beta$. This conjecture has been proven under a strong physical assumption that the time-ordered correlators factorize. As a consequence of that assumption, the magnitude of the normalized OTOC is always bounded by unity in the analytic regime on the complex time domain (see Refs. Maldacena et al. (2016); Tsuji et al. (2018) for detailed discussions). As another remarkable application of the OTOC-LE connection, we note that the factorization assumption can be removed for the Haar averaged OTOC, which is attributed to the fact that, in the corresponding complex time domain, the LE appears to be averaged over a thermal state with a positive effective temperature. We refer the interested reader to Appendix C for a detailed discussion.

Summary.— We have demonstrated the connection between two distinct areas of the dynamical quantum chaos, namely, the emerging field of the out-of-time-order correlators and the relatively more developed field of the Loschmidt echo. This relation not only allows a more general understanding of the universal properties of the OTOC but also provides new insights into both subjects. Two models were used to reveal the OTOC-LE relation in the scrambling regime and the intermediate exponential decay regime. Implications to the bound on the decay rate were also discussed. Future research could generalize the connection to higher-dimensional systems with spatial structures.

Appendix A Haar integral of subsystem operators

In this appendix we prove the formula in Eq. (4) in the main text.

For a single system, the Haar integral over a unitary operator $A$ gives

[TABLE]

Here, $dA\equiv d\mu(A)$ with $\mu$ the Haar measure on a unitary group.

This is a consequence of the defining property of the Haar measure, namely, the Haar measure is invariant under the transformation $A\rightarrow AU$, i.e., $\mu(A)=\mu(AU)$, where $U$ is an arbitrary unitary operator. This implies

[TABLE]

Thus, the Haar average of a operator $A$ is proportional to the identity operator. In finite dimension, its trace can be computed as

[TABLE]

The above integral is unique up to a constant multiplication factor; and the unitary groups on finite dimensional Hilbert spaces have finite measures. We can then choose $\int d\mu=1$ as a convention. The desired property in Eq. (13) is then deduced. We note that the identity $\int dA\ A^{\dagger}OA=\operatorname{tr}(O)\mathbb{I}$ holds in infinite dimensions as well. The formal mathematical treatment of the infinite dimensional case will be published elsewhere Bhattacharyya et al. (2019).

For an operator $O_{AB}$ on a composed system, the integral over operators on a subsystem can be performed in a similar way, by performing the partial trace:

[TABLE]

Appendix B Reduced dynamics for operator $B$

This appendix presents a more rigorous derivation of the reduced dynamical equation for the operator $B$ in Eq. (7) of the main text.

The total system Hamiltonian admits a general decomposition

[TABLE]

where $A$ denotes a small local subsystem $S_{A}$, $B$ denotes the compliment of $S_{A}$ to the total system. It is assumed that the total system is much larger than the local system $S_{A}$. $d_{A}$ is the dimension of $S_{A}$. The operators $\{V_{A}^{i}\}_{i}$ are chosen to be Hermitian and orthonormal, with respect to the (weighted) Hilbert-Schmidt inner product, i.e.,

[TABLE]

Here, $d_{A}$ is introduced as a convention to normalize the noise correlation functions (see below). If the operators $\{V^{i}_{A}\}_{i}$ are not normalized in this way, proper normalization factors will appear in the stochastic field as pre-factors. The operators $V_{B}^{i}$ on $S_{B}$ are also Hermitian, but their (Hilbert-Schmidt) norms are fixed and are equal to the norms of $H_{B}$. Thus, the parameter $\delta$ quantifies the relative strength of the couplings compared with $H_{B}$.

We are interested in strongly coupled systems, where the energy scales admits a hierarchy $\bar{H}_{A}\ll\bar{H}_{I}\ll\bar{H}_{B}$. For instance, in an $N$-particle system with all-to-all two-body interactions, when the subsystem $S_{A}$ refers to a single particle, the energy scales of $S_{A}$, $S_{B}$, and the coupling between them, are of the order of $1$, $N^{2}$ and $N$, respectively. Alternatively speaking, this condition means that the parameter $\delta\ll 1$; and energy scales of the subsystem $S_{A}$ is much smaller than that of the global dynamics.

Consider the reduced dynamics of an operator $B$ on the subsystem $S_{B}$ after the trace-out procedure. We have,

[TABLE]

This can be thought of as a decoherence process, i.e., the total system is prepared in an initial (unnormalized) product state $d_{A}\left(\frac{I_{A}}{d_{A}}\right)\otimes B$, where the subsystem $S_{B}$ is described by a “density matrix” $B$, and the subsystem $S_{A}$ is in a thermal state with infinite temperature up to normalization. The “quantum state” $B$ will become “mixed” with time evolution due to the presence of the couplings to subsystem $S_{A}$. When $\delta\ll 1$, the above evolution of $B(-t)$ can be expanded to the second order of $\delta$. This corresponds to the Born-Markov approximation, which leads the effective master equation for $B(-t)$ to a Lindblad form. It is known that in this case the effective master equation can be simulated with the evolution of $B$ under $H_{B}$ without coupling to other systems Budini et al. (1999), but subjects to a stochastic field

[TABLE]

with the correlations given by

[TABLE]

The approximation in the last step is due to the large energy hierarchy: the time scale of the dynamics of the subsystem $S_{A}$ is much larger than that of $B(t)$ under consideration. Alternatively, this can be thought of as taking the zeroth order the $H_{A}$. As a consequence, the noise field $l_{i}(t)$ can be taken as a random constant $\pm 1$ with equal probability. The reduced dynamics of the $B$ operator is then given by

[TABLE]

averaged over the stochastic field $\mathcal{F}$. Note that the pre-factor $d_{A}$ is present in the above equation due to the normalization of $\mathbb{I}_{A}$.

As the noise fields are random $\pm 1$, each realization of the stochastic field, denoted as $V_{\alpha}$, always appears as linear combination of $V^{i}_{B}$’s, whose coefficients take values $\pm 1$ randomly. The noisy evolution of $B(-t)$ is then given by

[TABLE]

where $N$ is the number of different realizations of the noise field.

Appendix C OTOC - LE connection

This appendix presents the derivation of Eq. (9) in the main text. It has been shown in the main text that after performing averaging over operator $A$, the OTOC reads

[TABLE]

Replacing the reduced dynamics of operator $B$ with its alternative expression Eq. (23) allows us to further evaluate the average over all unitary operators on subsystem $S_{B}$:

[TABLE]

This is Eq. (9a) in the main text.

The above equation can be generalized to finite temperatures, by distributing the operators at equal spacing around the thermal circle. Namely, let us define

[TABLE]

and evaluate the OTOC as

[TABLE]

By absorbing the inverse temperature into the time evolution operator, the OTOC can be evaluated in the same manner as in the infinite temperature case, but in the complex time domain, $t-i\beta/4$. The averaged OTOC is translated to an infinite temperature average of LE in complex time, which is further reduced to a finite temperature average of the regular LE in real time. More precisely (using Eq. (9b)),

[TABLE]

Here $O(V_{1},V_{2})$ is an operator resulting from the Baker-Campbell-Hausdorff expansion that contains terms involving the perturbation operators $V_{1}$ and $V_{2}$. This term is time-independent and only contributes to a vanishing correction of the thermal state. Note that the effective inverse temperature for the LE is $\beta/2$. The factor $Z^{2}_{\beta/2}/Z_{\beta}$ appears because of the particular regularization (27). As a consequence, the OTOC at $t=0$ has a non-unit value, even for unitary operators $A$ and $B$, in contrast to the regular thermal average scheme, where the OTOC at $t=0$ is unit for unitary operators. In our numerical simulations, as a convention, OTOCs are always normalized by removing this prefactor.

It is worth emphasizing that the Haar averaged OTOC is upper-bounded in the half strip $0<t$ and $-\frac{\beta}{4}\leq\tau\leq\frac{\beta}{4}$ for the complex time $t+i\tau$. The fact that the OTOC in this domain is bounded is crucial for the upper bound of the OTOC growth rate Maldacena et al. (2016); Tsuji et al. (2018). To see this fact, consider the averaged OTOC in complex time,

[TABLE]

which is a bounded quantity for the effective positive temperature $\beta/4\pm\tau$.

Appendix D OTOC and LE for the inverted harmonic oscillator

In this section, we give the details for the derivation of the OTOC and LE for the inverted harmonic oscillator.

D.1 Wavefunction of the IHO

The Hamiltonian of two coupled inverted harmonic oscillators (IHO) reads

[TABLE]

Here, for convenience, the frequencies are written as imaginary numbers. This Hamiltonian is reduced to the one of two uncoupled harmonic oscillators by the following transformation,

[TABLE]

where $\eta$ are $\xi$ are given by

[TABLE]

In the new coordinates, the effective masses and frequencies are given by

[TABLE]

The regularized thermal state evolved in the OTOC and LE is replaced with a pure product state (up to a normalization factor) $|\psi(x_{1})\rangle|\psi(x_{2})\rangle=e^{-x_{1}/m_{2}}e^{-x_{2}/m_{1}}$. The choice of this wavefunction is not physically essential for the OTOC and LE. It is however important from a practical point of view. This particular form of wavefunction is separable in the new coordinates $y_{1}$ and $y_{2}$, which makes solving the time-dependent wavefunction much easier.

For a given initial state $\psi(x)$ of a single IHO with mass $m$ and frequency $\omega$, the time-dependent wavefunction can be expressed by a path integral

[TABLE]

where the propagator is

[TABLE]

With these toolkits, the wavefunction $\psi(x_{1},x_{2},t)$ (which is involved in the OTOC evaluation; see the following section) can be expressed in the new coordinates using the single IHO propagator and then transfer back to the original coordinates.

D.2 Haar averaged OTOC for the IHO

In this section, we compute the Haar averaged OTOC for the IHO:

[TABLE]

The operators $A$ and $B$ to be averaged are unitary operators on the first and second IHO, respectively. Since a termal state of the IHO is not well defined, we replace the regularized thermal state $\rho_{\beta}$ with a pure state of the coupled IHO, i.e., $\rho(x_{1},x_{2};x_{1}^{\prime},x_{2}^{\prime})=\psi^{*}(x_{1},x_{2})\psi(x_{1}^{\prime},x_{2}^{\prime})=\rho^{\frac{1}{4}}(x_{1},x_{2};x_{1}^{\prime},x_{2}^{\prime})$. Using the property of the Haar measure described in Appendix A, the integral can be computed directly:

[TABLE]

Using the expression for the density operator, i.e.,

[TABLE]

the four terms in the above Haar integral can be evaluated separately as

[TABLE]

Besides the above Haar integral, the OTOC Eq. (36) involves a total trace, i.e.,

[TABLE]

After evaluating the trace and Haar integral, the OTOC finally reads

[TABLE]

Here, $\rho_{2}$ is the reduce density matrix of the second oscillator. The above equation is an exact result. In this particular setting, where the regularized thermal state is replaced with a pure state, the OTOC is precisely equal to the purity of one subsystem. The purity can be computed exactly as well, using the procedure described in the previous section.

D.3 Baker-Campbell-Hausdorff expansion of the LE

In this section we derive the second order expansion with respect to the perturbation strength $\delta$ for the Loschmidt echo of the type

[TABLE]

where $H$ is the Hamiltonian of a single IHO

[TABLE]

Here $\omega$ is a pure imaginary number. Treating $H_{0}\equiv H-\delta\hat{x}$ as the “unperturbed” Hamiltonian, the echo operator reads

[TABLE]

It can be expanded to the second order in $\delta$ using the Baker-Campbell-Hausdorff (BCH) formula Gorin et al. (2006):

[TABLE]

where $\Sigma(t)$ and $\Gamma(t)$ are determined by the Heisenberg evolution of the perturbation operator, i.e., $2\hat{x}(t)\equiv e^{iH_{0}t}2\hat{x}e^{-iH_{0}t}$,

[TABLE]

To solve the Heisenberg evolution of $2\hat{x}$ under Hamiltonian $H_{0}$, we perform a coordinate transformation $y=x+2\delta/(m\omega^{2})$. In the new coordinate the Hamiltonian reads

[TABLE]

This is a shifted harmonic oscillator Hamiltonian, for which the Heisenberg evolution of the position operator has been computed Hashimoto et al. (2017):

[TABLE]

which gives

[TABLE]

The last term in the above equation will be dropped since it contributes to higher than second-order in $\delta$ for the LE. $\Sigma$ and $\Gamma$ are thus given by

[TABLE]

In the following, we perform a Taylor expansion of the LE to second order in $\delta$. The second order expansion for the average of the echo operator (44) is given by

[TABLE]

Note that both $\Sigma(t)$ and $\Gamma(t)$ are real, there are only non-vanishing two second-order terms kept in the magnitude-square of the above quantity, which gives the final expression of the LE:

[TABLE]

As defined in the previous sections, the pure state for the OTOC average is chosen as a Gaussian form with zero mean. In this case only the first term in $\Sigma(t)$ has a non-zero contribution to the LE $M(t)$, i.e.,

[TABLE]

Appendix E Finite size study of the SYK model.

To further test our theory, we present numerical study of the fermionic version of the SYK model proposed in Fu and Sachdev (2016). The Hamiltonian of the SYK model reads

[TABLE]

where $J_{i,j;k,l}$ are complex Gaussian random couplings with zero mean obeying certain symmetries. $c_{i}$ and $c_{i}^{\dagger}$ are fermionic annihilation and creation operators at site $i$. We compute the OTOC for operators $c_{i}^{\dagger}+c_{i}$ on two distinct sites. Clear Gaussian decay with very weak temperature dependence has been observed (Fig. E.1), which agrees with previous numerical studies Fu and Sachdev (2016); Hosur et al. (2016) while contradicting the expected exponential decay with upper bound $2\pi/\beta$ Maldacena et al. (2016); Tsuji et al. (2018).

On the other hand, this observation fits into the theory of the present work: The OTOC decay rate is governed by the coupling strength between the target subsystem $S_{A}$, $S_{B}$ and the rest of the system. In a finite system, the coupling might be too large for the decay to be in the exponential regime. The decay rate extracted from our numerical simulation is close to the bandwidth (the width of the spectral density) of the SYK Hamiltonian, which indicates that the decay is in the Gaussian regime.

Since the relative strength of the coupling between the two subsystems $S_{A}$ and $S_{B}$ compared with the total Hamiltonian decreases with the system size, it is then expected that the OTOC decay would drop into the exponential regime in the large-$N$ limit. Due to the limited numerical capacity, we did not observe exponential decay up to 17 fermions, which is the largest system we are able to simulate numerically. However, we are able to observe exponential decay by manually adjusting the coupling strength. It is done by decreasing the couplings $J_{i,j;k,l}$ that involve the subsystem $S_{A}$ and $S_{B}$ by a factor $0<g<1$, while keeping the coupling in the rest of the system unchanged.

Figure E.2 shows the exponential OTOC decay for the deformed coupling strength. The decay rate also admits quadratic dependence on the coupling factor $g$, which satisfies the Fermi’s golden rule prediction for the exponential decay of the LE.