Out-of-time-order correlator in coupled harmonic oscillators

Exponential growth of thermal out-of-time-order correlator (OTOC) is an indicator of a possible gravity dual, and a simple toy quantum model showing the growth is being looked for. We consider a system of two harmonic oscillators coupled nonlinearly with each other, and numerically observe that the thermal OTOC grows exponentially in time. The system is well-known to be classically chaotic, and is a reduction of Yang-Mills-Higgs theory. The exponential growth is certified because the growth exponent (quantum Lyapunov exponent) of the thermal OTOC is well matched with the classical Lyapunov exponent, including their energy/temperature dependence. Even in the presence of the exponential growth in the OTOC, the energy level spacings are not sufficient to judge a Wigner distribution, hence the OTOC is a better indicator of quantum chaos.


Introduction
The out-of-time-order correlator (OTOC) [1] is considered as an indicator of quantum chaos. As it is a quantum analogue of the classical sensitivity against tiny changes in initial conditions, the exponential growth of the OTOC in its time evolution is regarded as the indicator of the chaos. The popular traditional indicator is the statistics of the energy level spacings (see for example [2]), and the interest is in whether the OTOC can be a better indicator of quantum chaos or not.
The OTOC attracted attention also because it can be an indicator of possible gravity dual, through the renowned AdS/CFT correspondence [3]. Gedanken experiments about shock waves in black hole geometries [4,5] have led to a maximum bound [6] of the quantum Lyapunov exponent, the exponent showing up in thermal OTOCs. The saturated value is interpreted as the red shift near the event horizon of a black hole with the Hawking temperature. The quantum Lyapunov exponent of the thermal OTOC can discriminate whether the system allows a gravitational picture. Theories with any gravity dual are strongly quantum, so we need a quantum analogue of the Lyapunov exponent as a discriminator, which is now provided by the thermal OTOC.
It was discovered [8] that the Sachdev-Ye-Kitaev (SYK) model [10,11], a quantum mechanical model of Majorana fermions, develops the exponential growth in its thermal OTOC and the quantum Lyapunov exponent saturates the bound. Stimulated by this discovery, research on OTOCs in various quantum systems have been carried out. Generic scheme for measuring the thermal OTOC in quantum mechanics was provided in [12] and it was shown that the thermal OTOC for a quantum stadium billiard, a typical chaotic system, does not grow exponentially. With adequate quantum states prepared, the expectation value of the OTOC was observed to grow exponentially for a kicked rotor [13], a stadium billiard [14], and the Dicke model [15], among few-body systems 1 . However, these examples are not with the thermal OTOCs. The thermal OTOC is indispensable for the identification with a gravity dual, hence a simple quantum system whose thermal OTOC grows exponentially is demanded.
In this paper, we study one of the simplest models: two harmonic oscillators coupled nonlinearly with each other, and show that the thermal OTOC of this model grows exponentially. Harmonic oscillators form a fundamental basis for any quantum field theory, which shows the importance of this simple system. In numerous examples, nonlinear couplings between harmonic oscillators provide classical chaos, needless to mention the popular example of the double pendulum. In our coupled harmonic oscillator (CHO) model, the chosen coupling x 2 y 2 with a coupling constant g 0 is significant by the following three reasons. First, the coupling g 0 x 2 y 2 is a part of the BFSS matrix theory [27] which is dual to a quantum gravity in the large N limit. As the coupling is a dimensional reduction of the Yang-Mills coupling, the model (1.1) serves as the starting point of a road to the AdS/CFT correspondence 2 . Second, this model (1.1) has attracted attention for long years [42][43][44][45] due to its classical chaos. It is one of the most popular chaos models, with known classical Lyapunov exponent and quantum level statistics [46,47]. Third, the model is obtained by a reduction from SU (2) Yang-Mills-Higgs theory which is a basis of the Standard Model of elementary particles 3 , so the model (1.1) is fundamentally related to the actual phenomena in the universe 4 .
We carefully demonstrate the numerical calculations of the thermal OTOC of the coupled harmonic oscillators (1.1), and see the exponential growth, by identifying the growth rate (named "quantum Lyapunov exponent") with the classical Lyapunov exponent of the system. This quantum-classical correspondence in the thermal OTOC, shown in the simple Hamiltonian system (1.1), may be extended to various other models with classical chaos, which will provide a novel arena for measuring quantum chaos.
This paper is organized as follows. In Sec. 2, we review the classical chaos of the coupled harmonic oscillator model (1.1). Sec. 3 is a review of [12] summarizing the numerical methods for calculating the thermal OTOC. Then in Sec. 4, we numerically evaluate the thermal OTOC of the system (1.1) and show that it grows exponentially in time. The quantum Lyapunov exponent measured there has the temperature dependence expected from the classical picture, as studied in ditail in Sec. 5. In Sec. 6 we show that at the energy scale of our concern the quantum energy level spacings are not sufficient to judge the quantum chaos, so the thermal OTOC is a better indicator of the quantum chaos. Sec. 7 is devoted to a summary and discussions. In App. A we estimate numerical errors in the evaluation of the OTOCs.

Review: Classical chaos in the coupled harmonic oscillators
In this section we introduce the coupled harmonic oscillator (CHO) model [42][43][44], which is the system we study throughout this paper, and make a review of its classical properties on chaos. In Sec. 2.1, we define the Hamiltonian of the model and describe the dimensional reduction. Then, in Sec. 2.2 we look at the Poincaré sections and the classical Lyapunov exponent, whose energy dependence is E 1/4 at high energy E [55].

The coupled harmonic oscillator model
The quantum mechanical model which we study is defined by the Hamiltonian (1.1), which we call CHO Hamiltonian: H = p 2 x + p 2 y + U (x, y) with the potential term Here ω and g 0 are constant parameters. Without losing its generality, we can take ω = 1 by a rescaling of the dynamical variables x(t) and y(t) and the Hamiltonian. In this paper, we take g 0 = 1/10, for our purpose of observing the chaotic and regular phases quantum mechanically 5 . This is a coupled harmonic oscillator system, and the coupling is given by the non-linear term g 0 x 2 y 2 .
To motivate readers, we here show that this Hamiltonian (1.1) can be obtained by a dimensional reduction of a four-dimensional SU(2) Yang-Mills-Higgs theory in a flat space-time 6 . The Lagrangian of the latter is given by The Levi-Civita symbol abc is a totally antisymmetric tensor with 123 = 1. The scalar field potential V (φ) in (2.2) is Although any quantum property of the Yang-Mills theory is quite intricate, after the dimensional reduction it becomes a tractable quantum model. First, we choose a gauge condition A a 0 = 0, and assume that the fields are spatially homogeneous; Here k, l denote spacial indices. In the vacuum which is spontaneously broken with µ 2 < 0, we turn on only A 1 1 ≡ x(t), A 2 2 ≡ y(t), then the Hamiltonian of the system reduces to where v ( = 0) corresponds to the vacuum expectation value of the original scalar field φ.
A simple rescaling brings this Hamiltonian to (1.1). It is known that, regardless of whether classical or quantum, this system (1.1) is in a regular phase at low energy and in a chaotic phase at high energy. This property is easily understood from the structure of the potential U (x, y). On the (x, y)-plane, in the region close to the origin, the potential is well approximated by the first term in U (x, y) in (2.1). Thus at the low energy, this system behaves as a two dimensional harmonic oscillator, and is in the regular phase. On the other hand, at the high energy, the energy contribution of the non-linear term becomes significant and this system turns into chaotic.
The following classical and quantum analyses have validated these phases. In the classical analysis [43,46,55], the Poincaré sections and Lyapunov exponents, which are useful to see if a classical system is chaotic or not, were examined. We review these quantities in the next subsection.
Quantum mechanically, the distribution of nearest-neighbor spacings of the energy eigenvalues can discriminate chaoticity of the system. Regular systems show a Poisson distribution, while chaotic systems show a Wigner one. In the quantum analysis [47], the distribution of the energy eigenvalues of the CHO model (1.1) was shown to be Wigner-like (Poisson-like) at high (low) energy.
In quantum systems, Out-of-Time-Order Correlator (OTOC) is another quantity to discriminate a quantum chaos. The quantitative indicator of quantum chaos is the quantum Lyapunov exponent, which measures the exponential growth of the OTOC. In this paper we focus on the OTOC of the CHO system (1.1), and evaluate the ability of it as an indicator of the quantum chaos, compared to the energy eigenvalue distribution method.

Energetic view of Poincaré sections and classical Lyapunov exponents
Drawing Poincaré sections is a popular method to find qualitatively a classical chaos. To illustrate concretely the classical chaos of the CHO Hamiltonian system (1.1), following [43,46], we plot in Fig.1 the Poincaré sections at several chosen values of the energy. In the figure, subfigure (a), (b), (c) and (d) correspond to E = 1, 2, 4 and 20 respectively. At a low energy, we see the orbits in the Poincaré section, so the system is in a regular phase. On the other hand, at high energy, orbits are destroyed and the Poincaré section is filled with scattered plots, meaning that the system is in a chaos phase. When one increases the energy from zero, the decay of the orbits starts to show up at E 2, and the whole Poincaré section is almost covered by the scattered plots at E 5. Hence the classical CHO Hamiltonian system has a transition between the regular phase and the chaos phase.
To be more precise, in our system (1.1) with the chosen parameters ω = 1 and g 0 = 1/10, the transition energy scale between the regular and the mixed phases is around E 2, and the one between the mixed and the complete chaos phase is around E 10. In the last section, we will compare our results of the quantum chaotic indicator OTOC with this classical phase structure.
The classical Lyapunov exponent λ cl measures how much a given classical system is chaotic. It characterizes the initial value sensitivity of the orbits in the phase space as an exponential growth ∆x(t) ∆x(0)e λ cl t , where ∆x(t) is the distance between two nearby orbits. Any positive λ cl means that the system is chaotic, and furthermore, the Lyapunov exponent can quantitatively evaluate the chaos. In Fig. 2 we show a set of numerically Figure 2. The classical Lyapunov exponents of our CHO Hamiltonian system (1.1) as a function of the energy of the initial condition. Note that we randomly chose the initial conditions for the numerical calculation of the Lyapunov exponents, and did not evaluate the maximum value of them for a given energy. This figure is for a rough understanding of the classical behavior of chaos of our system, and for the estimation of the value of the Lyapunov exponents. calculated 7 classical Lyapunov exponents of the CHO Hamiltonian system (1.1). At a low energy we find λ cl vanishes, which means that the system is regular. On the other hand, at a larger value of the energy the Lyapunov exponent is nonzero, thus chaos shows up. The measured Lyapunov exponent is an increasing function of energy, so the chaos is stronger for larger values of the energy. The energy dependence of the Lyapunov exponents is consistent with what we observed by using the Poincaré sections.
In fact, when the potential energy is just given by the interaction term, U (x, y) = g 0 x 2 y 2 , the energy dependence of the classical Lyapunov exponent can be analytically derived [55]: Ignoring the mass term (ω 2 /4)(x 2 + y 2 ) in the potential energy U (x, y) is equivalent to going to the high energy limit, so we can expect that at high energy our system (1.1) approximately follows the energy dependence (2.6). The relation (2.6) is derived as follows.
Ignoring the mass term in the Hamiltonian (1.1), we have the energy This system allows the following scaling transformation 8 : with an arbitrary positive constant λ. Since the Lyapunov exponent has the dimension of inverse time, from the last two relations in (2.8), it is clear that it should scale as (2.6) as a function of energy. 7 We used the mathematica LCE package formulated by M. Sandri. 8 The scaling law is valid even for the whole Yang-Mills theories and their supersymmetric generalizations.
In particular, the scaling law for the D0-brane matrix theory [27] was found in [53,54] and contributed to the early development in the study of M-theory.

Review: numerical computation of OTOC
In this section we summarize methods to numerically calculate OTOCs, based on [12]. First in Sec. 3.1, we define the OTOC which we focus on throughout this paper, and we verify that analytic expression of the OTOC for a harmonic oscillator is available, as its Heisenberg equation is exactly solvable. In general systems, however, it is difficult to solve the Heisenberg equation analytically, so we have to rely on a numerical method. We review the numerical methods for evaluating the OTOCs in generic quantum mechanical systems in Sec. 3.2.

Microcanonical OTOC and thermal OTOC
Let us take a quantum system which is described by a time-independent Hamiltonian H = H(x 1 , · · · , x N , p 1 , · · · , p N ), in the Heisenberg picture. In the following we write x ≡ x 1 , p ≡ p 1 , and denote the x operator in the Heisenberg representation at time t as x(t) ≡ e iHt x e −iHt . The OTOC of our concern is defined as Here · · · denotes the thermal average with the inverse temperature β ≡ 1/T , Assuming discrete energy eigenvalues, we have the following expression by using energy eigenvectors, where H|n = E n |n with a non-negative integer n. We call C T (t) a "thermal OTOC" and c n (t) a "microcanonical OTOC." In a quantum chaotic system, it is expected that the thermal OTOC grows exponentially as C T (t) e 2λq(T )t , where λ q (T ), which we call "quantum Lyapunov exponent," is positive. This exponent is a quantum analogue of the classical Lyapunov exponent described in the previous section.
As a concrete example, we consider a harmonic oscillator in one dimension, In this case it is possible to solve the Heisenberg equation analytically, and the solution is Substituting these solutions and using the canonical commutation relation, we get A higher dimensional harmonic oscillator has the same expression since the degrees of freedom in the other dimensions are decoupled from the sector of (x 1 , p 1 ). The OTOCs of the harmonic oscillator are oscillatory, thus do not show any exponential growth. This is consistent with the fact that this system is solvable and regular. In addition, they depend on neither the mode label n nor the temperature T .

Numerical method for OTOC
In general quantum systems it is difficult to solve their Heisenberg equation, so we rely on numerical methods to evaluate the OTOCs. Here, we rewrite the OTOCs (3.3) for our later purposes to performing the numerical calculations. The completeness condition n |n n| = 1 brings the microcanonical OTOC c n (t) in (3.3) to the following form: Substituting x(t) = e iHt xe −iHt to this expression of b nm (t) and using the completeness condition again, we obtain Here we have defined x nm ≡ n|x|m , p nm ≡ n|p|m , and E nm ≡ E n − E m . When the Hamiltonian is of the form (which includes the case of the CHO Hamiltonian (1.1)) one can easily verify [H, x] = −2ip. Then evaluating the matrix elements of this relation with n| and |m , we find p nm = i 2 E nm x nm . By using this, we derive a useful expression With these, OTOCs can be numerically evaluated, following the concrete steps: 1. Solve the Schrödinger equation of the given system to obtain the energy eigenvalues and the wave functions.
2. Compute x nm = n|x|m with numerical integration.
3. Substitute the result of 2 to (3.10) to calculate b nm (t).
4. Evaluate the microcanonical OTOC c n (t) by substituting the result of 3 to (3.7).
In these numerical evaluations, approximations by introducing a finite cut-off to the infinite sums are necessary 9 . In the next section, we shall use this strategy to numerically evaluate the OTOCs in the CHO system.

Preparations: microcanonical OTOC
As a preparation for evaluating the thermal OTOC from the next subsection, here we describe our results of the microcanonical OTOC obtained through the numerical strategy 1,2,3 and 4 given in the previous section. First, we solve the time-independent Schrödinger equation 10 and obtain the energy eigenvalues E n and the wave functions ψ n (x, y). See Fig. 3 for the obtained distribution of the energy levels. By using them, we evaluate the microcanonical OTOC. Fig. 4 is our numerical result. The truncation errors in this evaluation will be dealt with in Sec. A.
The resultant microcanonical OTOC for low modes oscillates almost periodically in time, and any exponential growth of it cannot be found. This behavior is similar to the microcanonical OTOC of a harmonic oscillator (3.6). On the other hand, for higher modes, the behavior of the microcanonical OTOC deviates from that of the harmonic oscillator. This is expected from the structure of the potential U (x, y) of (1.1). The potential U (x, y) is well approximated by a harmonic oscillator around the origin of the (x, y) plane, so the wave functions localized in the vicinity of the origin behave like that of the harmonic oscillator, which happens at low energy. On the contrary, for the higher modes, the the wave functions spread, and the second term of the potential U (x, y) contributes dominantly. This explains the deviation of the microcanonical OTOC from the harmonic oscillator (3.6).

Reading quantum Lyapunov exponent from thermal OTOC
By substituting the numerical results of the microcanonical OTOC c n (t) to the expression (3.3), we compute the thermal OTOC C T (t). The numerical results are shown in Fig. 5. The obtained thermal OTOC shows two apparent characteristics: at the early time stage 10 For solving it, we use Mathematica package NDEigensystem. it grows, and at the late time stage it saturates to a constant value. In Sec. 4.2.1 we study that the latter part reproduces the expected property in any quantum chaotic system. Then in Sec. 4.2.3 we show that the growth of the thermal OTOC at early time is exponential.

Reproducing late-time behavior of OTOC
It is known that in generic quantum chaotic systems, the correlation between x(t) and p is gradually lost in the time evolution due to the chaoticity. After the Ehrenfest time, the value of the thermal OTOC (3.1) is expected to approach asymptotically in time [6] We shall use this first for checking our numerical evaluation of the thermal OTOC. Let us look at our numerical results of the thermal OTOC, Fig. 5. At low temperature, the thermal OTOC does not grow so much, as in the case of the low mode microcanonical OTOC. This is because, in (3.3) at low temperature, the contribution from high modes is suppressed due to the Boltzmann factor e −βEn . On the other hand, at higher temperature, the Boltzmann weight does not suppress the high modes, then more number of modes contribute to the thermal OTOC. Indeed, our numerical results in Fig. 5 show that at late times the time evolution of the value of the OTOCs is flattened, and look like they approach some constants.
To find out quantitatively the asymptotic values, we evaluate the value of the OTOC for each temperature by averaging its values at t = 20, 22, 24, 26, 28, 30. In Fig. 6 we plot those asymptotic values evaluated from our numerical results of the thermal OTOC, as a function of the temperature. For a comparison, we also plot the temperature dependence  of the expected expression (4.2). Their trends are quite similar to each other. Thus we conclude that the saturated parts of the thermal OTOC indicate that the CHO Hamiltonian system is quantum chaotic in the temperature region we consider 11 .

Very early stage of OTOC
Next we focus on the early time part of the thermal OTOC in Fig. 5. At the very early stage, the thermal OTOC decreases, and the temperature dependence is not observed. Then after a while the temperature dependence shows up. To extract the quantum Lyapunov exponent due to the quantum chaos later, we shall reveal the physical origin of this behavior at the very early stage.
At sufficiently early times, the contribution of the nonlinear term 1 10 x 2 y 2 of CHO Hamiltonian (1.1) can be treated as a perturbation; with g 0 = 1/10 and ω = 1. For that reason, the very early time behavior should be identical to that of the harmonic oscillator. There exists a time scale at which the nonlinear term starts to contribute significantly. This time scale can be estimated in the following way. For small t, the Heisenberg operator x(t) = e iHt xe −iHt in the OTOC (3.1) can be approximated by . Here x 0 (t) is the Heisenberg operator evolved by the free Hamiltonian H 0 , which is given in (3.5). We adopt a principle that when the expectation value of g 0 x 2 y 2 t for the ground state |0 of H 0 grows roughly to O(10%), the perturbation expansion given above is broken. Evaluating this with the harmonic oscillator ground state, we find 0|g 0 x 2 y 2 t|0 = g 0 t , (4.4) 11 The slight difference from (4.2) could come from our truncation of energy levels, and also from the finite time effect.
therefore we expect that, since our coupling constant is g 0 = 1/10, t 1 is the time when the nonlinear term starts to contribute to the value of the thermal OTOC. This estimation is consistent with the numerical results shown in Fig. 5.

Reading quantum Lyapunov exponent
We are ready to read the quantum Lyapunov exponent from the thermal OTOC. Our numerical results in Fig. 5 show that the time evolution behavior of the OTOC is complicated, but we can use the observations in Sec. 4.2.1 and Sec. 4.2.2 to extract the time regions for which the OTOC is expected to evolve exponentially in time. As seen in Sec. 2.2, the classical CHO system is subject to a phase transition from the regular to the chaotic phase, when going from lower to higher energy scales. Equating the energy roughly to the temperature, the "classical phase transition temperature" is about T 3. So, it is expected that the thermal OTOC can grow exponentially for T ≥ 3. Furthermore, in Sec. 4.2.1 and Sec. 4.2.2, we observe that for t ≤ 2 the nonlinear effect is not strong enough, while for t ≥ 6 (which is the expected Ehrenfest time) the system saturates to the asymptotic value of the OTOC. From these observations, we make a numerical fitting of the data of our thermal OTOC for T ≥ 3 in the time domain 12 2 ≤ t ≤ 6, by an exponential function a(T )e 2λq(T )t with the temperature-dependent quantum Lyapunov exponent λ q (T ). We plot our numerical result for the quantum Lyapunov exponent for various values of the temperature, T = 3, 3.5, 4, 4.5, 5, in Fig. 7.
Due to the following two observations about the exponents, we claim that the thermal OTOCs of the CHO system grow exponentially.
(1) The order of magnitude of λ q equals that of the classical λ cl .
From the first place, the OTOC C T (t) is defined as a quantum analogue of the classical Poisson bracket {x(t), p} 2 e 2λ cl t , therefore the behavior C T (t) e 2λ cl t is expected, as long as the system is still close to the classical time evolution, until Ehrenfest time. Thus it is also natural to expect that the order of quantum Lyapunov exponent λ q and of classical Lyapunov exponent λ cl is same. Our classical Lyapunov exponent are shown in Fig. 2 and the values 13 are O(0.1) ∼ O(0.2), which is consistent with the 12 To judge whether the evolution is linear or exponential, we need the time period which is not much shorter than the inverse of the Lyapunov exponent. Our time period is δt = 6 − 2 = 4, while we will see that our quantum Lyapunov exponent of the thermal OTOC gives 1/λ ∼ 4, which is almost equal to the time period. Also, note that for the case of the quantum stadium billiard [12], the growth of the thermal OTOC was not judged to be exponential, because even when one takes a deformation to the regular limit (by bringing the billiard to a circle), the growth remains. In our case, if one takes the deformation to the regular limit by bringing g0 → 0, the system reduces to the decoupled harmonic oscillators so our growth in the thermal OTOC disappears. From this observation, it is plausible that the growth in our thermal OTOC is exponential and is due to the chaos. We would like to thank K. Murata for his valuable comment on this. 13 Note that the plots in Fig. 2 are just for some fixed initial conditions, and we haven't searched the largest Lyapunov exponent. So the maximal values should be a little bit larger than the values in Fig. 2. Furthermore, the quantum Lyapunov exponent is for the thermal average, which includes not only the mode at E = T but also the other energy modes. Thus equating these two exponents need a caution. Here we just check whether the order of magnitude for these two exponents is shared commonly or not. quantum values in Fig. 7.
(2) The temperature dependence of λ q is similar to the energy dependence of λ cl .
We postulate that the temperature dependence of the quantum Lyapunov exponent λ q is λ q (T ) = bT c with some constants b and c, and fit our numerical results with this function, see Fig. 7. The fitting function satisfies λ q (T = 0) = 0, since at zero temperature there should not be the exponential growth 14 . As a result, we obtain 15 On the other hand, as shown in (2.6) in Sec. 2.2, the energy dependence of the classical Lyapunov exponent is given as λ cl ∝ E 1/4 at sufficiently high energy. Regarding the typical energy at finite temperature is E T , we expect that the quantum Lyapunov exponent λ q (T ) is approximately proportional to T 1/4 . Our result (4.6) is consistent with that 16 .
The coincidence with the classical behavior is remarkable. Although in the CHO system which has only two degrees of freedom the Ehrenfest time is expected to be rather short, we find a time domain in which the classical exponential growth can be detected in thermal 14 By the way, at T → 0, only ground state of the microcanonical OTOC contributes the thermal OTOC because of Boltzmann factor e −βEn . The wavefunction of ground state is localized in the region close to the origin on (x, y) plane, and is almost the same as a harmonic oscillator's one. Therefore 5) and the thermal OTOC is not expected to grow exponentially. It is natural to assume λ = 0 at T = 0. 15 The error is a systematic error: the exponent c depends on how many points we adopt when we fit them. 16 The assumption E T is further verified in details, in Sec. 5. OTOC. In other words, when the classical behavior such as the energy dependence of the Lyapunov exponent is known, one can identify where the quantum thermal OTOC behaves classically, and estimate the Ehrenfest time.
One of our major results is that the temperature dependence of the quantum Lyapunov exponent of the CHO system is found as (4.6). The temperature dependence is important as it indicates how close the system is to a black hole system [6]. For the exponent of generic OTOCs 17 there exists a bound 2λ q (T ) ≤ 2πT in the large N limit. The measured quantum Lyapunov exponent λ q (T ) in our system is found to be consistent with this bound.

Temperature and energy dependence of Lyapunov exponent
In the identification of the quantum Lyapunov exponent in Sec. 4.2.3, we used an argument that the temperature dependence of the quantum Lyapunov exponent λ q (T ) is expected to be λ q (T ) ∝ T 1/4 due to the energy dependence (2.6) of the classical Lyapunov exponent λ cl (E). In fact, this could be a too rough assumption -as discussed in [14], the thermal average can wash off quantum chaotic information encoded in the microcanonical OTOCs. Here in this section, we elaborate on this concern in detail. We evaluate the thermal average and show by a numerical calculation that the assumption we made in Sec. 4

.2.3 is reasonable.
Since the issue is the relation between the microcanonical OTOC and the thermal OTOC, and that between classical and quantum correlations, first we rewrite the thermal OTOC (3.3) as The definition of the OTOC here is different from that of [6] as for the ordering of the Boltzmann weight operators. Here we naively compare the bound with ours, by assuming that the bound is valid also for our OTOC of the CHO system (3.1).  Here, we have defined where N (E) counts the degeneracy of the states at energy E. We set c E (t) = 0 if there exists no energy level at E. Next, we approximate the discrete energy level distribution of the CHO Hamiltonian system by a continuous function. We assume that c E (t) corresponds to the classical Poisson bracket {x(t), p} 2 P , which means In Sec.2.2 we have seen that at high energy the classical Lyapunov exponent λ cl (E) shows the energy dependence (2.6). We use it here as a simple assumption that the relation (2.6) holds for all energy scales of our concern, and set In what follows, we numerically calculate (5.1) and investigate its temperature dependence, to show that the assumption made in Sec. 4.2.3 is a reasonable one. In the first place, we need to find the density of states ρ(E) for calculating (5.1). For that purpose, we define W (E) which is the number of states below a given value E of the energy. We numerically solve the Schrödinger equation for the CHO Hamiltonian system up to E 40. to energy E, we find the density of states ρ(E) as Using this, the thermal OTOC (5.1) is 18 (5.7) Next, we make the E integration of (5.7), for which we need the normalization value a of the classical Lyapunov exponent. We adopt the value a = 0.20 obtained by a numerical evaluation of the classical Lyapunov exponent at E = 40 (which is high enough so that the system is totally chaotic as seen in Fig. 1) 19 .
Considering the time domain t ≤ t Max = 6 which is our expected Ehrenfest time observed in Fig. 5, and defining a rescaled time variable τ ≡ at (and τ Max = at Max = 1.20 accordingly), we perform the integral (5.7) by discretizing 20 possible values of τ and T .
The numerical result C T (τ ) of (5.7) at each fixed T is fit by an exponential function αe Lτ , see Fig. 9. We find 21 the exponent L as a function of T , see Fig. 10. Then finally we numerically fit this set of L's by a power function of T , γT c with constants c and γ, as 18 Note that (5.6) is valid only for E 40 where the Schrödinger equation is solved. With (5.6), we can deal with the thermal OTOCs at a temperature value at which the states with E 40 dominate. In Sec. 4, we numerically calculated the thermal OTOC for T ≤ 5. When T = 5, the Boltzmann factor e −βE for the energy E = 40 is e −8 0.00034, which is small enough. Thus, when T ≤ 5, the contribution of the energy levels with E ≥ 40 is well suppressed. In the numerical calculation of (5.1), we consider T ≤ 5 and ignore E ≥ Etrunc = 40. 19 The classical Lyapunov exponent numerically depends on the initial conditions in the phase space. We will evaluate the effects of this deviation from a = 0.20 later. 20 We discretize the variable τ by units of 0.01: τ = 0.01m (m = 0, 1, · · · , mMax = 120), and the temperature T by units of 0.1: T = 0.1n (n = 0, 1, · · · , 50). 21 Note that this L is different from what we call the quantum Lyapunov exponent by the time-rescaling factor a. However, since our target is just the power of T in the quantum Lyapunov exponent, we don't need to worry about the overall rescaling factor. Then we numerically find that the exponent c takes slightly different values depending on τ Max , see Fig. 11. The resultant c is distributed in a rather small range 0.26 < c < 0.28. So we conclude that the result (5.8) is trustable against the errors in evaluation of the classical Lyapunov exponent.

Comparison to energy level statistics
So far, we have investigated the OTOC of the CHO model (1.1). In Sec. 4 we found the exponential growth of the thermal OTOC, which indicates the quantum chaos, for the temperature range 3 ≤ T ≤ 5. On the other hand, as we mentioned in Sec. 2.1, the distribution of the nearest-neighbor spacings of the energy eigenvalues is often used to discriminate chaoticity of the system. Regular systems show a Poisson distribution, while chaotic systems show a Wigner one. In the quantum analysis [47], the distribution of the energy eigenvalues of the CHO model was shown to be Wigner-like (Poisson-like) at high (low) energy.
We are interested in whether our thermal OTOC is a better indicator of the quantum chaos compared to the energy level statistics. The analysis made in [47] used the CHO model with parameters different from ours. So, for the comparison let us find a relation between the parameters of ours and those of [47]. Our Schrödinger equation is given in (4.1) with ω = 1 and g 0 = 1/10, while that of [47] is The numerical calculations in [47] were made with the choice k = 1/200. When the nonlinear term is dominant and the system is totally chaotic, we may ignore the harmonic potential terms 22 . In that approximation, for (6.1) to be identical to (4.1), we need to make the rescaling x → ax and bẼ = E, with a = (5/2) 1/6 and b = 2 2/3 5 1/3 . So we find the relation between the energy E of our system and the energyẼ of [47], In the analysis of [47], the energy level spacings show the Wigner distribution for the energỹ E ≥ 85, while show the Poisson one forẼ ≤ 85. Translated to our energy, this threshold energy scale is E 230. Now, our energy levels used for the evaluation of the thermal OTOC is E 32. This means that the detection of the chaos by means of the OTOCs is much more effective compared to the statistics of the energy level spacings.
To make sure that the argument above is correct, here we explicitly analyze the energy level spacings of our system (1.1). The statistical analysis needs to take care of the symmetry of the system, which is the point group C 4v in our case. We consider wave functions which are completely symmetric under the C 4v group. First, to compare our OTOC analyses, we consider the energy eigenstates which we used for the OTOC analyses: E 32 (see Appendix A for the truncation scheme). Fig. 12 Left is the obtained histogram for the nearest-neighbor spacings of energy eigenstates with the eigenvalue E )/(35 − 1). It is obvious that we lack the number of states sufficient to judge whether the histogram is subject to the Wigner distribution or not.
To see the distribution, we need more states. If we use 100 states which are symmetric under C 4v (whose highest energy goes up to E ∼ 60), we obtain a histogram shown in Therefore, we conclude that by the level statistics it is almost impossible to see the Wigner distribution only with the states which we used for the evaluation of the OTOCs.
Quantitatively, we can evaluate whether a given distribution is Wigner-like or not by ther-parameter [59,60], which is defined as r n = min(s n , s n−1 ) max(s n , s n−1 ) , s n = E (sym) If we did not ignore the harmonic potential term, the exact comparison between our (4.1) and (6.1) could not be made. So, here, we adopt that approximation. We expect that even with this approximation the order of estimate is still valid.
E or T ≈ Figure 13. Summary table of the energy (temperature) dependence of the classical/quantum chaos of the CHO system (1.1). The thermal OTOC can detect the quantum chaos at lower energy (temperature) compared with the level statistics. The quantum regular/chaos phase structure through the OTOC is found to be consistent with the classical phase structure, which indicates that the thermal OTOC properly reflects the property of the classical system.
By definition,r n is a real number between 0 and 1. The distribution is characterized by the average ofr-parameter, which we denote as r . For example, r is around 0.386 for a Poisson distribution. For a Wigner distribution, r takes a value in 0.531 ∼ 0.674. For the 35 energy eigenvalues in Fig. 3 Left, we find r 0.478 , (6.4) which is different from that of a Wigner distribution. For the 100 energy eigenvalues in Fig. 3 Right, we find r 0.526 , (6.5) which is closer to the region for the Wigner distribution. Thus, the distribution of nearest-neighbor spacings of the energy eigenvalues which we used for the calculation of the OTOC is not sufficient to judge a Wigner distribution. In other words, we cannot discriminate quantum chaoticity of the CHO model by the level statistics in the energy scale. In spite of this, we have observed the exponential growth of the thermal OTOC in Sec. 4. We can conclude that the OTOC is more sensitive to quantum chaoticity than the level statistics, and is a better indicator of quantum chaos.

Conclusion and discussion
In this paper, we numerically calculated the microcanonical OTOC and the thermal OTOC of the nonlinearly coupled harmonic oscillator (1.1), for various temperature and energy range. The nonlinear coupling is a reduction of SU(2) Yang-Mills theory which is expected to have a gravity dual in the large N and strong coupling limit. We confirmed the exponential growth of the thermal OTOC at higher temperature, and could read the quantum Lyapunov exponent. The Lyapunov exponent measured has a temperature dependence λ q ∝ T c with c = 0.26 ∼ 0.31.
We have two reasons for the identification of the quantum Lyapunov exponent. First, the order of magnitude of the calculated exponent (Fig. 7) coincided with that of the classical Lyapunov exponent, λ cl (Fig. 2), as shown in Sec. 4. Second, the temperature dependence of the calculated exponent (4.6) reproduces the dependence (5.8) expected from the energy dependence of the classical Lyapunov exponent, λ cl , (2.6), as shown in Sec. 5.
In Sec. 6 we showed that the level statistics of the energy eigenstates that are effective for the calculation of the OTOC does not show a Wigner distribution. This means that the level statistics fails in discriminating quantum chaos of the coupled harmonic oscillator model at the energy scale of our concern. On the other hand, we found the exponential growth of the thermal OTOC. Therefore the thermal OTOC is a better indicator of quantum chaos than the level statistics, for the nonlinearly coupled harmonic oscillator system. Fig. 13 summarizes our quantum results and their relationship with the classical regular/chaos phase structure. The horizontal axis is the temperature (or energy) of the quantum and classical systems. In the classical coupled harmonic oscillator system described in Sec. 2.2, at low energy, we see the orbits in the Poincaré section, so the system is in a regular phase, while at high energy, orbits decay and the Poincaré section is filled with scattered points, meaning that the system is in a chaos phase. In the quantum system, the thermal OTOC is more sensitive to the quantum chaos than the level statistics, and is able to discriminate the quantum chaos at lower energy (temperature). Moreover, the "phase transition" temperature between the OTOCs oscillatory phase and exponential growth phase is observed to close to the energy at which the classical regular/chaos phase transition occurs. Thus the thermal OTOC properly reflects the chaotic property of the classical system.
It is encouraging that even the very simple coupled harmonic oscillators exhibit the exponential growth of the thermal OTOC. Inclusion of more oscillators will make the Ehrenfest time later, which will help reading the quantum Lyapunov exponent. For the system to get closer to the black hole system, according to the generic AdS/CFT dictionary, we need to take the large N limit and the strong coupling limit. For the latter we need to go to lower temperature because the coupling constant has a positive mass dimension, and thus need also to suppress the harmonic oscillator frequency ω to have the chaotic phase at lower energy scales. Achieving these limits will enables us to find the quantum Lyapunov exponent saturating the bound 2λ q = 2πT , the black hole behavior. The emergence of this truly quantum regime would be quite intriguing for understanding quantum gravity.
Also in the introduction we emphasized the relation to the holographic principle, that is why we considered the thermal OTOCs. Another motivation for using the thermal OTOC is that it captures a global feature of the phase space. Classical instability could be easily captured by localized wave functions with which the expectation value of the OTOC is taken. However, if one localizes the wave function near the top of the potential hill, then even when the whole system is not chaotic, the OTOC can grow exponentially (see for example [61]). The local growth does not immediately mean chaos, and in fact, this kind of exponential growth can occur even in an inverted harmonic oscillator [62] or in a doublewell potential [63] in one dimension. The thermal OTOCs, on the other hand, use energy eigenstates which are spread over the whole phase space, thus would not see these local instabilities. Since ergodicity is expected for classically chaotic systems, even widely spread wave functions for the thermal OTOCs can detect the chaos due to the global structure of the phase space.
One of the importance of our system (1.1) is its genericity. Harmonic oscillators are everywhere, as a basis of quantum field theories and quantum information theories. They are also easier to be implemented as experimental apparatus, and at finite temperature the thermal OTOC is appropriate to serve as a better diagnosis for finding quantum chaos. We expect that some harmonic oscillator couplings other than our (1.1) will produce quantum chaos similarly, with different scaling properties. Studying the thermal OTOCs and scrambling timescales in systems given by some limits of coupled harmonic oscillators will open the way to explicitly examine quantum gravity nature of matters. Based on the considerations above, for computational reason, we choose N trunc = 253, and take 150 for the number of the microcanonical OTOCs c n (t) for calculating the thermal OTOC. This N trunc = 253 corresponds to the state with its energy E 32.