Out-of-time-order correlators in quantum mechanics

The out-of-time-order correlator (OTOC) is considered as a measure of quantum chaos. We formulate how to calculate the OTOC for quantum mechanics with a general Hamiltonian. We demonstrate explicit calculations of OTOCs for a harmonic oscillator, a particle in a one-dimensional box, a circle billiard and stadium billiards. For the first two cases, OTOCs are periodic in time because of their commensurable energy spectra. For the circle and stadium billiards, they are not recursive but saturate to constant values which are linear in temperature. Although the stadium billiard is a typical example of the classical chaos, an expected exponential growth of the OTOC is not found. We also discuss the classical limit of the OTOC. Analysis of a time evolution of a wavepacket in a box shows that the OTOC can deviate from its classical value at a time much earlier than the Ehrenfest time.

a vertex correction of a current for a superconductor [1], was recently turned to be considered as a measure of the magnitude of quantum chaos. A naive argument for the relation between the OTOC and chaos is as follows [2]. Consider position and momentum operators, x(t) and p(t), in a quantum system. We can define an OTOC as C T = − [x(t), p(0)] 2 . Taking a naive semiclassical limit, we would be able to replace the commutator [x(t), p(0)] by the Poisson bracket i {x(t), p(0)} PB = i δx(t)/δx(0). For a classically chaotic system with a Lyapunov exponent λ, we have δx(t)/δx(0) ∼ e λt because of the sensitivity to initial conditions. Thus, the OTOC should grow as ∼ 2 e 2λt and we can read off the quantum Lyapunov exponent λ from it. The quantization of a classically chaotic system may provide a positive quantum Lyapunov exponent of the OTOC. A possible distinction from the classical chaotic system is that the OTOC does not grow eternally but saturates at the Ehrenfest time t E . The Ehrenfest time is defined by the time scale beyond which the wave function spreads over the whole system. It is roughly characterized as a boundary between a particle-like behavior and a wave-like behavior of the wave function.
In recent years, the OTOC has been regarded as an important observable in the context of AdS/CFT correspondence [3] or quantum gravity. A maximum bound of the quantum Lyapunov exponent was proposed as λ ≤ 2πk B T / [2]. The bound was originally suggested in the context of quantum information around black hole horizons [4,5] (see also Refs. [6][7][8][9][10]). The Lyapunov bound is saturated by the Sachdev-Ye-Kitaev (SYK) model [11,12]: A quantum mechanics of Majorana fermions with infinitely long range disorder interactions. Saturation of the quantum Lyapunov bound indicates that the SYK model describes a quantum black hole through the AdS/CFT correspondence.
Since the original calculation of the quantum Lyapunov exponent by Kitaev, there appeared subsequent study for generalizing the SYK model [13][14][15]. However, we are still missing explicit examples of OTOCs. Do typical chaotic systems show exponential growth in OTOCs? Can we find any qualitative difference between integrable and chaotic systems through OTOCs? To answer these problems, we study the OTOC of single particle quantum mechanics. First we formulate how to calculate the OTOC for generic quantum mechanics. In particular, by the reason described above, we choose W = x and V = p to measure a possible indication of quantum chaos. Based on the formalism, we examine OTOCs of some popular quantum systems: (i)a harmonic oscillator, (ii)a particle in a one-dimensional box, (iii)a circle billiard (a particle in a circle-shaped infinite well), and (iv)a stadium billiard. Former three are known as integrable systems. The stadium billiard [16][17][18][19][20], on the other hand, is one of the most popular and well-studied Hamiltonian chaotic systems.  T represents the temperature of the systems. We find a clear distinction between the two: The OTOC for the particle in a box periodically comes back to its initial value (= 1), while that for stadium billiard saturates to a constant value. The asymptotic value grows with the temperature. In the inset of the right panel, we show the early time evolution of the OTOC. We find no clear exponential growth of the OTOC.
Among our main results, we show two numerical results in Fig. 1, which shows typical behavior of the OTOCs. The left panel is a numerical evaluation of the OTOC for a particle in a 1D box, and the right is that for a stadium billiard. In the figure, we took the unit of = k B = 2m = 1 where m is the mass of the particle. We also set (Length of the box) = 1 or (Area of the billiard) = 1. We summarize our findings in this paper below: 1. The OTOCs grow at early times. However, at least for T 400, they do not show apparent exponential growth even for the stadium billiards.
2. The OTOC of a particle in a box is periodic in time because of its commensurable energy spectrum, while that of the stadium billiards saturates to a constant value.
3. The OTOCs grow with temperature except for the harmonic oscillator. The high temperature limit does not reproduce the classical value in general. For high temperature, asymptotic or maximum values of OTOCs grow linearly in temperature as C T ∼ mT × (typical system size) 2 .
4. An analysis of a time evolution of a wave packet in a 1D box shows that the OTOC deviates from its classical value at the time scale parametrically earlier than the Ehrenfest time t E .

Organization of this paper
We start in Sec. 1 with an introduction and a summary of our results obtained in this paper. We formulate how to calculate the OTOCs in a general quantum mechanical system in Sec. 2. Then in Sec. 3 we evaluate the OTOCs for integrable examples such as the harmonic oscillator, a particle in a 1D box, and a particle in a circular billiard. In Sec. 4, after reviewing the classical chaos of the stadium billiards, we present our numerical results of the OTOCs for the quantum stadium billiards. In Sec. 5, we study time evolution of a wave packet in a 1D box, to find a deviation of the OTOC from its classical value at rather early times. Sec. 6 is devoted for discussions, with a relation to quantum fidelity and Loschmidt echo. Our appendices include description on our numerical truncation errors and detailed formulas for the analytic calculation of the OTOC for a wave packet.

Units in this paper
Throughout this paper, we work with the unit of = k B = 2m = 1, where m is the mass of a particle. When we consider a particle in the 1D box and the stadium billiard, we also set (length of the box) ≡ L = 1 and (area of the billiard) ≡ A = 1, respectively. For the billiard, one can easily restore dimensional parameters notifying For the particle in the box, A is replaced by L 2 .

Out-of-time-order correlators in quantum mechanics
In this section, we propose a formalism to compute the OTOC for general quantum mechanics with time-independent Hamiltonian: H = H(x 1 , · · · , x n , p 1 , · · · , p n ). We consider the out-of-time-order correlator (OTOC) defined by Here we define β = 1/T with the temperature of the system T . We denoted x = x 1 and p = p 1 for notational simplicity. Hereafter, we will omit the argument of Heisenberg operators for t = 0: O ≡ O(0). Taking energy eigenstates as the basis of the Hilbert space, we can rewrite the OTOC as where H|n = E n |n . We will refer the OTOC for a fixed energy eigenstate, c n (t), as a microcanonical OTOC. On the other hand, we will refer C T (t) as a thermal OTOC. Once we compute microcanonical OTOCs, we can obtain the thermal OTOC by taking their thermal average. 2 Let us rewrite the microcanonical OTOC using matrix elements of x and p for numerical calculations. Using the completeness relation 1 = m |m m|, we rewrite the microcanonical OTOC as . Substituting x(t) = e iHt xe −iHt and inserting the completeness relation again, we obtain Although OTOCs for the harmonic oscillator have been already studied in Ref. [21], we compute them again using the formalism in the previous section. The energy spectrum and matrix elements of x is given by where n, m = 0, 1, 2, · · · . Substituting above expressions into Eq.(2.7), we have b nm (t) = δ nm cos ωt . They are periodic functions whose periodicity is ∆t = π/ω. They do not depend on energy level n or temperature T . We will find that this is a special property only for the harmonic oscillator amongst our examples. As in Ref. [21], one can also get the same result using the explicit expression of the Heisenberg operators: From the explicit solution, we have [x(t), p(0)] = i cos ωt and OTOCs are given as Eq. (3.4). This method is not useful for other cases since it is difficult (or impossible) to obtain explicit expressions of Heisenberg operators for a general Hamiltonian.

Particle in a box
One of the other integrable examples is a particle in a 1D box. The Hamiltonian for the one-dimensional case is Eigenfunctions and eigenvalues are given by where n = 1, 2, · · · . The matrix elements of x are written as Although we know exact expressions of x and energy eigenvalues, it would be impossible to carry out the summation in Eq.(2.7) analytically. So, we evaluate b nm (t) Microcanonical out-of-time-order correlators for the particle in a 1D box (L = 1). and c n (t) numerically with truncation n, m ≤ N trunc = 100 and compute OTOCs. In Figs.2 and 1(a), we show microcanonical and thermal OTOCs, respectively. Note that the energy spectrum for the particle in a box is commensurable: E n is proportional to the integer n 2 . By using Eq. (3.8), one can show that the all E kl appearing in rhs of Eq. (2.7) become π 2 × odd integer and thus, OTOCs have periodicity ∆t = 1/π. For large n, microcanonical OTOCs become large and tend to oscillate since high frequency modes become relevant in Eq.(2.7). High frequency oscillations seem to be suppressed in the thermal OTOC. The thermal OTOC also tends to be large at high temperature. We have found that the maximum of the thermal OTOC increases linearly as a function of T : max C T ≃ 0.1672 × 2mT L 2 (mT L 2 ≫ 1) where the size of the box L and the mass of the particle m are restored. We have also checked that the time average of the OTOC is given bȳ whereC T = lim τ →∞ τ 0 dtC T (t)/τ . For the particle in a 2D square box, V (x, y) = 0 (0 < x < 1, 0 < y < 1), ∞(else), we can obtain the same result for the thermal OTOC. Eigenstates in the 2D box are completely separable as |n x |n y , so the operator [x(t), p(0)] 2 does not operate to |n y . Therefore, the existence of the y-direction is completely irrelevant for calculating the thermal OTOC.

Circle billiard
As a non-trivial 2D example, we consider a circle billiard: In this case, the x-and y-directions are not separable unlike the 2D square box. It is known that classical dynamics of a particle in the circle billiard is integrable. (We will see that in section.4.1.) We fix the radius of the circle as R = 1/ √ π so that the area of the billiard becomes unity. Taking polar coordinates x = r cos θ and y = r sin θ, we obtain exact expressions for eigenvalues and eigenfunctions as where k ∈ Z and l ∈ {1, 2, · · · }. J k is the Bessel function of the first kind and ρ kl represents its l-th root, i.e. J k (ρ kl ) = 0. The normalization factor is given by The energy spectrum for the circle billiard is not commensurable. It is only asymptotically commensurable: It tends to be commensurable for high energy because of ρ kl ≃ (k/2 + l)π for l ≫ 1. Although eigenvalues and functions are labeled by two integers k and l, we relabel them by a single integer n in ascending order of E kl and denote them as (E n , ψ n ). The matrix elements of x can be obtained from We can carry out the integration of θ analytically. We perform the numerical integration along r-direction and obtain the matrix elements. Substituting the matrix elements and energy eigenvalues into Eq.(2.7) and using Eq.(2.3), we obtain OTOCs. Fig.3 shows the microcanonical and thermal OTOCs for the circle billiard. The microcanonical OTOCs seem to be non-periodic and tend to be larger for a larger energy level n. We can also find "dips" in the microcanonical OTOCs: For example, for n = 40 and 100, they become small (c n ∼ O(1)) around at t ≃ 0.8 and t ≃ 1.4, respectively. In the thermal OTOCs, we can also find similar dips around at t ≃ 0.8. We will see that, for the stadium billiard, the dip does not appear in OTOCs. The dips in OTOCs would originate from the asymptotically commensurable property of the spectrum and be reflecting the integrability of the systems.

Non-integrable example: Stadium billiards 4.1 Classical mechanics of stadium billiards
As a typical example of the non-integrable (chaotic) system, we consider a stadium billiard [16][17][18][19][20]: The domain Ω is shown in Fig.4(a). We denote radii of semicircles as R and the length of straight lines as 2a. Let us revisit the classical dynamics of the particle in  the billiard. Inside the stadium, the particle moves freely with a constant velocity. At the boundary of the stadium, the particle is reflected elastically. Fig.4(a) also shows a typical trajectory of the particle in the stadium. We can find the chaotic behavior.
One of the most characteristic behavior in chaotic systems is the sensitivity to initial conditions: A tiny deviation of the initial condition causes a significant difference in the future. The Lyapunov exponent is a useful quantity to measure the strength of the sensitivity to initial conditions. Denoting the phase space variable as X(t), we consider its linear perturbation: X(t) → X(t) + δ(t). If X(t) is a chaotic solution, because of the sensitivity to initial conditions, the perturbation expands exponentially as δ(t) ∼ e λt . The growing rate λ is called Lyapunov exponent. A positive Lyapunov exponent is the signal of chaos.
In Fig.4(b), we show the Lyapunov exponent as a function of the deformation parameter a/R. 3 Here, we took the unit of v = A = 1, where v is the velocity of the particle and A = πR 2 + 4aR is the area of the stadium. From the dimensional analysis, we can easily restore v and A by replacing λ → √ Aλ/v. The Lyapunov exponent is zero at the circle limit a/R = 0. Hence, the classical circle billiard is integrable. For positive a/R, λ increases quickly and reaches maximum value around at a/R ∼ 1.3. The rough estimation of the Lyapunov exponent is where v and A are restored. In case of the dynamical billiard, the Lyapunov exponent is proportional to the velocity v, apparently. (The frequency of collisions is proportional to the velocity.) In the squeezed limit a/R → ∞, the particle does not have any chance to hit the semicircles of the stadium. Thus, λ also approaches zero  The maximum Lyapunov exponent vs deformation parameter a/R for fixed area of the stadium and velocity of the particle, A = v = 1.
in this limit. This result is consistent with earlier calculations of Lyapunov exponents in Refs. [20,22,23].

Quantum mechanics of stadium billiards
As the quantum version of the dynamical billiard, we consider the time-independent Schrödinger equation [−∆ + V stad (x, y)]ψ n = E n ψ n [24]. To determine the eigenvalues and eigenfunctions, we used the Mathematica standard package, NDEigensystem, which solves eigenvalue problems of linear differential operators using the finite element method. In Fig.5, we plot eigenvalues of the quantum stadium billiard with a/R = 1. The energy spectrum is roughly linear in the energy level n. Fitting the spectrum, we have E n ≃ 13n (a/R = 1). This approximate formula is useful for rough estimation of the energy spectrum. In Fig.6, we show eigenfunctions of the quantum billiard with a/R = 1. It is known that the Eherenfest time, at which the wave function spreads over the whole system, becomes quite small for the chaotic system [25][26][27][28]. To illustrate it, we consider the macroscopic billiard: m = 1 kg, A = 1 m 2 and v = 1 m/s. For these parameters, from Eq.(4.2), the Lyapunov exponent is estimated as λ ∼ 1 Hz. Such billiard seems sufficiently classical, but actually, there is tiny uncertainty in its position and momentum. Let us take the uncertainty as ∆x ∼ 10 −17 m and ∆p ∼ 10 −17 kg m/s so that the uncertainty principle is saturated: ∆x∆p ∼ . In the chaotic system, the wave packet of the particle would exponentially spread as ℓ(t) ∼ ∆xe λt . When the size of the wave packet becomes the same order as the system size, L ∼ 1 m, a quantum interference effect becomes significant. The Eherenfest time is estimated as t E ∼ λ −1 ln(L/∆x) ∼ 40 s. So even if we start from the extremely localized wave packet, just after one minute, the system becomes completely quantum. This behavior is different from what we find in nature. The problem was that we assumed that the system is isolated from the environment. Once we take into account the weak interaction between the system and environment, decoherence is caused and "the decoherence suppresses the quantum suppression of the chaos" [25]. For instance, the emergence of the classical chaos due to the decoherence is discussed by considering the continuous quantum measurement [29].
In this paper, we consider the isolated quantum billiard. Even for high temperature or high energy, after the Eherenfest time, the quantum effects will be important and classical approximation will breakdown. We will revisit this point in section 6.

Out-of-time-order correlators
From the eigenfunctions, we obtain matrix elements of x as x nm = Ω dxdy ψ n xψ m . Substituting x nm and E n into Eq.(2.7) and using Eq.(2.3), we compute the micro- canonical OTOCs as functions of t for each energy level n. In Fig.7, we show the microcanonical OTOCs for the stadium billiard with a/R = 1. For n = 1, 2, OTOCs look similar to those for the particle in a box. (See Fig.2.) This is because typical scales of the wave functions for small n are of the same size as that of the system. So, wave functions do not "feel" the curvature of semicircles of the stadium. For higher n, however, OTOCs become less recursive than that for the circle billiard and oscillate around constant values at late time. Taking the thermal average of the microcanonical OTOCs, we compute the thermal OTOC. In Fig,1(b), we show the thermal OTOC for the stadium billiard. For low temperatures, the lower n mode dominates the thermal OTOC and it looks similar to the microcanonical OTOC for n = 1. For high temperature, the thermal OTOC increases quickly as a function of t and approaches a constant value at late time. The magnitude of the oscillation around the constant value is small compared to the OTOC of the circle billiard.
In particular, we do not observe dips found in the circle billiard. We have done same calculations for a/R = 0.2i (i = 1, 2, · · · , 10) and found qualitatively similar behavior.
Can we find an exponential growth in the OTOCs? Fig.8 shows an early time evolution of thermal OTOCs for stadium and circle billiards. The OTOC for the stadium billiard does not show a clear exponential growth. At a very early time t 0.01, one may be able to argue that there is an exponential region. However, to find the exponential growth e λt , we need much longer time than 1/λ. (Otherwise, we cannot distinguish the exponential and linear functions.) Moreover, a similarlylooking behavior can be found even for the circle billiard. There is no qualitative difference in early time OTOCs between the stadium and the circle billiards. Our results indicate that, at least for T 400, we cannot distinguish integrable and chaotic systems from the early time evolution of the thermal OTOCs. In Ref. [2], it was proposed that the Lyapunov exponent λ defined by C T (t) ∼ e 2λt satisfies a bound λ ≤ 2πT . The thermal OTOC of the stadium billiard does not show the exponential growth and, in that sense, it trivially satisfies the bound. We can observe that thermal OTOCs approach constant values at late times. What determines the asymptotic value? A naive expectation is that the OTOC saturates when it becomes the "system size". Since the OTOC has the dimension of 2 , the asymptotic value of the OTOC would be given by C T ∼ P 2 sys L 2 sys where P sys and L sys are the typical momentum and size of the system. For a thermal system with the temperature T , the typical momentum would be P sys ∼ √ mT where m is the mass of the particle. Therefore, our expectation is We can numerically confirm this relation for the stadium billiard as follows. We evaluate the asymptotic values of thermal OTOCs from C T (t = ∞) ≃ t 2 t 2 dtC T (t)/(t 2 −t 1 ). We took t 1 = 5 and t 2 = 10 in actual calculations. Fig.9(a) shows C T (t = ∞) as functions of T for several choices of the deformation parameter a/R of the stadium shape. Our numerical results clearly show that C T (t = ∞) linearly depends on T and its slope depends on a/R. In Fig.9(b), we plot the slope C T (t = ∞)/T as function of a/R. It is also given by a linear function of a/R. Fitting the plot, we obtain where the area of the billiard A is restored. Substituting A = πR 2 + 4aR, we can rewrite above expression as C T (∞) ≃ 0.68(a + 0.94R)(a + 0.79R)mT . Since the system size of the stadium is given by L sys ∼ 2(a + R), this is consistent with the naive prediction from the dimensional analysis (4.3). In the limit of R → 0, the system reduces to the particle in a 1D box with L = 2a. For R → 0, we obtain C T (t = ∞) ≃ 0.0858 × 2m(2a) 2 T . This is also consistent with the time average of the OTOC for the particle in the box (3.9). First we discuss that the classical statistics does not reproduce the high temperature limit of the OTOCs in general. Using the example of the particle in the 1D box, we can easily show that classical statistics is not so useful for estimation of the OTOC.
In the high temperature limit, a naive expectation is that the the thermal average in the OTOC can be replaced by the integral in the 2D phase space as where Z cl = dxdp/(2π) e −βH and { , } PB is the Poisson bracket. For the particle in the box, the classical solution is explicitly written as before the bounce at a boundary. After the bounce, the momentum is reflected as p(t) → −p(t). We consider the infinitesimal deviation of the initial position fixing the momentum as (x(0), p(0)) → (x(0) + δx(0), p(0)). By the time evolution, the particle will bounce at boundaries. Then, the deviation of the position change its signature but the absolute value is constant: δx(t) = (−1) n δx(0) after n-th bounce. Therefore, we have Substituting this into Eq.(5.1), we obtain This classical value is apparently different from the quantum result of the OTOC at a high temperature shown in Fig.1(a). We can also estimate the classical OTOC for the stadium billiard. From the sensitivity to initial conditions, we have {x(t), p(0)} PB ∼ e λt . From Eq.(4.2), the Lyapunov exponent is λ ∼ v ∼ p(0) for A = 1.
Although, for fine-tuned initial conditions, the particle motion can integrable, their measure would be zero. For t ≫ β, we can replace ∞ 0 dp by ∞ −∞ dp and we have This has unusual dependence in t and is again apparently different from the quantum calculations in Fig.1(b).
In the case of the harmonic oscillator, on the other hand, the classical solutions, x(t) and p(t), are completely identical to Eq.(3.5). Therefore, classical statistics gives the same result as the quantum calculation: C

Out-of-time-order correlator for a wavepacket
Why does quantum statistics not approach classical statistics? To answer the question, we consider a simpler setup: OTOC for a wavepacket in a 1D box. We will show that the OTOC deviates from its classical value at a time much earlier than the Ehrenfest time.
The wavefunction of the wavepacket is given by We consider the well localized wavepacket in real and momentum spaces: Expanding this wavepacket by the energy eigenstates (3.7), we obtain Using the wavepacket, we consider expectation values of commutator [x(t), p(0)] and its square as Here, b nm (t) has been defined in Eqs.(2.4) and (2.7). We know the analytic expression for the matrix element of x (3.8) and the energy spectrum (3.7) for the particle in a box. We perform the summation numerically. We set parameters in the wavepacket as k 0 = 3000π, σ = 1/(30π) and x 0 = 1/2. For |n − 3000| ≤ 120, we use the expression in Eq.(5.9) as α n . For |n − 3000| > 120, we simply set α n = 0. Fig.10 shows the time evolution of the wavepacket |φ(t, x)| 2 . At the early time, the wavepacket is well localized in the real space and shifts with a constant velocity v = k 0 /m = 2k 0 . The wavepacket is getting spread as time increases. The width of the wavepacket σ t spreads as 2) in the appendix.) When the width of the wavepacket is the same order as the system size, a quantum interference effect becomes important. Hence the classical (particle) interpretation is no longer valid. This time scale is called Ehrenfest time. For the particle in a box, the Ehrenfest time t E is estimated from σ t | t=t E ∼ (system size) ∼ 1 and we have t E ∼ σ .
around at the (2N + 1)-th bounce. Here, ℓ(t) = 2k 0 t + x 0 − (2N + 1) is the difference between the center of the wavepacket and the right boundary x = 1. Before the Ehrenfest time σ t ≪ 1, the 2-point OTOC is approximated by a step function. Therefore, for a well localized wavepacket in a box, we have a quantum-classical correspondence: On the other hand, for the 4-point OTOC c φ , we can observe the spiky profile at the time of the bounce: k 0 t = 0.25 + 0.5n (n = 1, 2, 3, · · · ). Except for the spiky points, it is well approximated by the classical prediction. By an analytic calculation in appendix.B, for x 0 = 1.2, we obtain around at the bounce at the boundary. Spikes in right panel of Fig.11 are gaussians whose widths are given by σ t . We focus on the time just on the bounce: ℓ(t) = 0. Then, from the inequality of arithmetic and geometric means, the spike term in above equation becomes 42.0 σk 2 0 + 6.38 The spike term grows to be the same order as the classical value c classical φ = 1 by This is sufficiently earlier than the Eherenfest time because of t s /t E ∼ 1/ √ k 0 σ ≪ 1. Therefore, for the wavepacket in the box, we would be able to say We need a shorter time scale t ≪ t s to see the quantum-classical correspondence in the 4-point OTOC.

Discussion
The OTOC of the stadium billiard does not show the exponential growth. We also found that classical statistics for the OTOC does not coincide with the quantum calculation even for the particle in a 1D box. For the discussion on the disagreement, we have to be careful about the classical limit of OTOCs because it is expected that the classical behavior shows up only at times earlier than the Ehrenfest time t E . Let us estimate the Ehrenfest time t E for the thermal system of the particle in a box. At a temperature T , the typical energy of the particle should be E ∼ T . So, the typical momentum of the particle is estimated Although there is no a priori choice for the typical size of a particle σ, it should satisfy T −1/2 ≪ σ ≪ 1 from the well-localized-condition (5.8).
We take the thermal de Broglie length as the typical size of the particle σ ∼ T −1/2 since this gives the smallest Ehrenfest time. Then, from Eq.(5.13), the Ehrenfest time is estimated as t E ∼ T −1/2 . For T ∼ 100, we have t E ≃ 0.1. In Fig.1(a), even if we focus on the time scale of t ≪ t E , the OTOC disagrees with the classical value C cl (t) = 1. We can argue the stadium billiard in a similar manner. For the chaotic system with a Lyapunov exponent λ, the width of a wavepacket would spread exponentially as σ t ∼ σe λt . The Ehrenfest time, which is estimated from σ t | t=t E ∼ 1, is given by t E ∼ λ −1 ln(σ −1 ). For a thermal system, the typical velocity is given by v = k 0 /m ∼ √ T . Then, from Eq.(4.2), the Lyapunov exponent is λ ∼ √ T . Choosing the typical size of the particle as σ ∼ T −1/2 again, we can estimate the Ehrenfest time as t E ∼ T −1/2 ln T . For T = 400, we have t E ∼ 0.3. In our numerical result of the OTOC given in Fig.8(a), we cannot find an exponential growth for the time region t ≪ t E .
Why do the OTOCs deviate from their classical value at a high temperature and t ≪ t E ? In section.5.2, we found the other time scale t s , at which the quantumclassical correspondence of the 4-point OTOC is violated, for a wavepacket in a box. We showed that t s is sufficiently smaller than the Ehrenfest time t E . Although we do not have any physical interpretation of t s at the moment, the existence of the time scale t s would be an origin of the distinction between quantum and classical mechanics as for the OTOCs. The time scale t s might stem from the interference effect at the bounce (Fig. 11). In fact, such a small time scale does not show up in Ref. [30], in which the system without the boundary was considered.
It has been also known for quantum fidelity or Loschmidt echo (see [31] for a review) that generally the time region for reproducing the classical Lyapunov behavior is quite limited. The Loschmidt echo measures how identical the state is to   Fig. 12, generic OTOCs are interpreted as a generalization of the Loschmidt echo. 4 Therefore, it is natural that the semiclassical limit of the OTOC of the billiard does not reproduce the classical Lyapunov behavior, as in the case of the fidelity.
From asymptotic values of OTOCs for the stadium billiards, we found the empirical relation for the typical magnitude of the thermal OTOC: C T ∼ mT × (system size) 2 . This result indicates that the magnitude of the OTOC does not relate to the magnitude of chaos. In fact, while the classical Lyapunov exponent has a maximum value around at a/R = 1.3 as in Fig.4(b), the magnitude of the OTOC is just given by a linear function in a/R for fixed A as in Eq. (4.4).
By a naive argument in Sec. 1, the OTOC can be related to the classical Lyapunov exponent via the replacement of the commutator by a Poisson bracket. However, in our analyses of the quantum stadium billiards, we do not find the exponential growth of the OTOC. Is there single particle quantum mechanics which shows clear exponential growth in the thermal OTOC? Can the quantum Lyapunov exponent saturate the bound provided in Ref. [2]? To answer these questions, we need further study of OTOCs of classically chaotic systems.
The OTOC for the Sachdev-Ye-Kitaev (SYK) model grows exponentially [11,12]. What was essential for the exponential growth? There are two significant difference between the SYK and our examples. (1)The OTOC in the SYK model has been calculated in a large N limit while our examples concern single particle quantum mechanics. For the large N theory, we can divide the system into two parts, A and B. The part B can be regarded as the "environment" by integrating out the degree of freedom in B. The interaction between the environment and part A would cause the decoherence [25][26][27][28]. It follows that the system would be classical-like and show the exponential growth in the OTOC. Indeed, the emergence of the decoherence by taking the partial trace for the environment have been shown [39,40]. The coupled systems, each of them classically shows the chaotic behavior, has also been investigated and shown to have decoherence effect [41]. (2)The SYK model has the random coupling. The random coupling is known to enhances the decoherence [42]. Adding to that, it is known that the wave function is localized in space when the system has random potential (Anderson localization). The localization of the wave function would be regarded as emergence of particle nature and thus the exponential growth of the OTOC might be expected since the presence of the wave nature in our model prevents us to observe the exponential growth. However, the localization is also known to have negative effect for classicalization since the diffusion is suppressed due to the localization in the momentum space [43]. Thus the effect of the randomness on OTOC is still unclear. n < 100, we found better convergence than n = 100. For n > 100, microcanonical OTOCs does not contribute to the thermal OTOC so much because of the suppression factor exp(−E n /T ). (In this paper, we consider T ≤ 400 for the stadium billiard. For n = 100, the energy eigenvalue is E 100 ≃ 1300 and its contribution is suppressed by exp(−E n /T ) ≃ 0.04.) Based on the analysis in this section, we chose N trunc = 400 for most of calculations of stadium billiards.

B.1 Propagation of a wavepacket in a box
We consider dynamics of a wavepacket in a 1D box: V (x) = 0 (0 < x < 1), ∞ (else). The initial gaussian wave function is given by Eq.(5.7). We consider the well localized wavepacket satisfying Eq.(5.8). We also assume that the center of the wavepacket is separated from boundaries: Then, we do not have to mind tiny non-zero values of the wave function at boundaries. For the free particle V (x) = 0, time evolution of the wavepacket is given by where α(t) = 1 + it/σ 2 and U free (t) is the time evolution operator of the free particle: U free (t) = e −ip 2 t . The absolute square of the wavepacket is given by a gaussian: where σ t is defined in Eq. (5.12). The width of the wavepacket spreads as a function of t. The center of the wavepacket is moving with a constant velocity v = 2k 0 . In case of the particle in a box, the dynamics of the wavepacket is given by the "folding operation" of the free wavepacket as . One can check that this satisfies Schrödinger equation i∂ t φ = −∂ 2 x φ and boundary conditions φ(t, x = 0, 1) = 0. Hereafter, we only consider much earlier time than the Ehrenfest time, For following calculations, it is convenient to introduce the "folding operator" F by Using the folding operator, the time evolution operator for the particle in a box is written as Also the Hermite conjugate of the time evolution operator is written as U † (t) = U(−t) = F U free (−t) = F U † free (t). One can easily check following formulae of the folding operator: where n ∈ Z.

B.2 Operation of x(t) and p(0) to the wavepacket
We consider around (2N + 1)-th bounce of the wavepacket at boundaries: The center of the free wavepacket x center = 2k 0 t+ x 0 is in 2N < x center < 2N + 2. We also assume that the wavepacket does not overlap with the left boundary: x center − 2N ≫ σ and 2N + 2 − x center ≫ σ. Then, dynamical solution can be approximated by where we define For the computation of the OTOC, we consider the operation of x(t) and p(0) to the wavepacket. For the initial gaussian wavepacket, we can rewrite the operation of the momentum operator as The operator A commutes with x and p since it does not contain x and ∂ x .
We consider the operation of the Heisenberg position operator x(t) − 1 to the wavepacket (We consider x(t) − 1 instead of x(t) for the simplicity of the following calculations.): At the last equality, we used Eq.(B.10). By the similar way as the momentum operator, the operation of (x − 1) to free wavepacket can be written by using k 0derivative as The operator B commutes with x and p. The introduced variable ℓ(t) represents the coordinate difference between the center of the free wavepacket x = x center and the position of (2N + 1)-th bounce x = 2N + 1. One can check that the introduced operators A and B satisfy the "canonical commutation relation": Using the operator B, we obtain We need to calculate the inverse time evolution of φ free (t, x ± )h(x). The strategy is same as the previous subsection: We consider the inverse time evolution by the free Hamiltonian U † free and apply the folding operator F . The propagator of the free particle is given by Using the propagator, we have At the last equality, we extended the lower bound of the integration to −∞ because we assumed that the wavepacket is not around the left boundary. Completing the square in the exponent of the integrand, we can rewrite above expression as (B.21) Note that we can rewrite ( 20). Using the error function 5 , we have Taking the folding operation, we obtain inverse time evolution of the wave function in a box as At the second equality, we used the definition of x ± (B.11) and the formula of folding operator (B.8). At the last equality, we used the other formula (B.9). So, we obtain where we used erf(−z) = −erf(z) and defined

B.3 2-point out-of-time-order correlator
We can easily obtain analytic expression for the 2-point OTOC. The 2-point OTOC is given by At the third equality, we replaced momentum operator by A. We already know the wave function |φ(t) as in Eq.(B.10). Therefore, the 2-point OTOC is written as We neglected the cross terms such as φ * (t, x + )φ(t, x − ) since they oscillate very quickly as ∼ e ±2ik 0 x and canceled out by the integration. Substituting the explicit expression of φ fee (t, x) (B.2) and introducing x ′ = x − 1, we obtain the 2-point OTOC as At the bounce, the 2-point OTOC changes the signature. Its time scale is given by ∆t ∼ σ t /k 0 . This is consistent with the numerical calculation in Fig.11(a).
In the curly bracket of Ψ L , −iα(t) − i/2 is negligible. We can see that as follows.

B.5 4-point Out-of-time-order correlator
The 4-point OTOC is given by (B.49) The cross terms of right and left movers are negligible in the integration. Recall that At the second equality, we used the fact that Φ(y) is suppressed by 1/q = 2σ/(y − x 0 ) ≪ 1 outside the domain of the box. At the third equality, we used Eq.(B.25).