Regularization dependence of the OTOC. Which Lyapunov spectrum is the physical one?

We study the contour dependence of the out-of-time-ordered correlation function (OTOC) both in weakly coupled field theory and in the Sachdev-Ye-Kitaev (SYK) model. We show that its value, including its Lyapunov spectrum, depends sensitively on the shape of the complex time contour. By choosing different contours one can violate the Maldacena-Shenker-Stanford-bound on chaos. Our result puts into question which of the Lyapunov exponents computed from the exponential growth of the OTOC reflects the actual physical dynamics of the system. We argue that, in a weakly coupled $\Phi^4$ theory, a kinetic theory argument indicates that the symmetric configuration of the time contour, namely the one for which the bound on chaos has been proven, has a proper interpretation in terms of dynamical chaos. This is supported by the SYK model with $q/2$-body interactions, where for $q\gg1$, and at the conformal point, the contour dependence disappears due to kinematical reasons, and its value is also that of the symmetric contour configuration. Finally, we point out that a relation between these OTOCs and a quantity which may be measured experimentally -- the Loschmidt echo -- also suggests a symmetric contour configuration, with the subtlety that the inverse periodicity in Euclidean time is half the physical temperature. In this interpretation the chaos bound reads $\lambda \leq \frac{2\pi}{\beta} = \pi T_{\text{physical}}$.

It has long been known that chaos, understood as the exponential sensitivity of the dynamics to initial conditions, does not have an immediate equivalent in the quantum dynamics governed by the Schrödinger equation. In quantum systems one needs to define quantum chaos in a more indirect way. One way to do so, is to measure the correlation between an operator W (t) and some earlier perturbation V (0) and compare this with the correlation where the perturbation V (0) is performed after operator W (t) is inserted: Choosing W (t) = q(t) and V (0) = p(0) this commutator formally equals [W (t), V (0)] = i ∂q(t) ∂q(0) and in that sense the above measures the sensitivity to initial conditions. The commutator is evaluated between two wave-functions, however. For a generic |ψ initial and |ψ final , this is a complex amplitude that also depends on the details of both. An obvious step is to sum over final states, which converts this to an expectation value To also isolate the dynamics driven by V (0) and W (t) as much from the details of the initial state, one can average over a suitable ensemble. A physically natural choice is the thermal one This commutator square C(t; β) or, equivalently, this out-of-time ordered correlation function (OTOC) has been of much interest as a diagnostic of chaotic behaviour in many-body systems [1][2][3]. Specifically, if this OTOC has a regime where it exhibits an exponential time dependence: C(t) ∼ e λt this behaviour has been proposed to be a signature of chaos, with λ being the quantum Lyapunov exponent. 1 Moreover, this quantum Lyapunov exponent has been conjectured to be bounded from above λ ≤ 2πk B T / [3].
In practice most computations do not compute C(t) as defined above. Rather one "smears" the thermal distribution between the two commutators [3,4] C(t; β) regulated ≡ Tr ρ  ( 1.4) Mathematically, this has the advantage of being manifestly Hermitian (see e.g. [4]).
The physical intuition is that in a QFT this correlation function naively suffers from a short-distance divergences caused by the insertion of two operators at the same time. As chaos is in principle a long-time characteristic, the claim is that the information about chaos, and in particular the Lyapunov exponents λ, do not depend on this regularization [3,4]. We will show that this intuition is incorrect, as was also pointed out earlier in [5] for the specific case of 2D fermions with quenched disorder. By explicit computation we will show that in the two-parameter family of "regularized" OTOCs 1 Note that the Lyupanov exponent defined this way is in fact twice the chaos exponent one would surmise from the choice W (t) = q(t), V (0) = p(0) with q(t) ∼ e λ chaos t q(0), i.e. λ = 2λ chaos .
the Lyapunov exponents are independent of σ but do depend on α. In fact, in the SYK model, for any α = 1/2 these Lyapunov exponents will violate the MSS-bound, though we stress that the bound was only derived for the specific case α = 1/2 [3]. Our computation shows that this regularization dependence is an IR-effect and has nothing to do with short-distance singularities. The more appropriate comparison for the regularization dependence of the OTOC is the proof in Schwinger-Keldysh theory that physical correlation functions are independent on the choice of contour. In Schwinger-Keldysh theory, there is a diagrammatic proof that physical Green's functions involving operator insertions either on only forward or only backward branches are independent of the contour due to energy conservation; this can be found in e.g. [6,7]. The OTOC, however, is a correlation on a doubled Schwinger-Keldysh contour [8] and the two-body Green's functions involved in the commutator square involve operators inserted on both forward and backward branches. The arguments of [6,7] do not generalize to prove that the correlation functions that appear in C(t; β) (α,σ) must be independent on the contour. Our explicit computation shows that they indeed are not. This is not a pedantic point, as also pointed out by [5]. OTOCs are now being measured either in numerical or actual physical experiments. Often one massages the regulator to be the most convenient for the set-up. For instance, Das et. al. [9] use the canonical thermal OTOC C(t; β) (0,0) in a numerical study, whereas a cold atom experiment measures a Loschmidt echo [10], which can be related to C(t; β) ( 1 2 ,0) . As the theoretical prediction for these two correlation functions is different due to the regulator dependence, these two experimental results cannot be compared to each other.
Given the regularization dependence that we and [5] observe, the immediate question arises: which is the proper regularization that measures quantum chaos. As the previous paragraph shows, to some extent this is in the eye of the beholder. One can devise experimental set-ups that measure either. Nevertheless, we will argue that the OTOC that most closely reflects physical microscopic chaos is the symmetrized one C(t; β) ( 1 2 ,0) used originally for hermiticity reasons. Our argument rests on two facts: 1. In weakly coupled field theories the computation of any of the OTOCs C(t; β) (α,σ) can be cast in the form of a kinetic equation [11]. This kinetic equation reveals most closely the physical process one is actually computing. In terms of the kinetic equation, only the symmetrized OTOC with α = 1/2 can be understood as a microscopic unbiased "collision"-counter. Such unbiased collision counters have long been successfully proposed as tracking microscopic classical chaos [12,13]. This is explained in section 3.2.
2. In strongly coupled Sachdev-Ye-Kitaev (SYK) models, the Lyapunov exponent also depends on the contour if computed straightforwardly. This therefore directly affects the Lyapunov spectrum near the emergent conformal regime at strong coupling. However, in SYK models with q/2-body interactions there exists a large-q approximation where contour dependence disappears due to kinematical reasons. We will explain this in section 4. Computing the IR Lyapunov exponent in the OTOC in this large q-approximation, the result qualitatively agrees with the OTOC computed with α = 1/2 and finite q, and should be exactly equal to that in the limit βJ 1.
We conclude by showing the symmetric OTOC C(t; β) ( 1 2 ,0) regulated this way has a natural interpretation as a Loschmidt echo, rather than an expectation value in a thermal ensemble as in the introductory thought experiment. This has as subtle physical consequence that the physical temperature is set by twice the inverse periodicity in Euclidean time. In this interpretation the MSS bound reads

A two-parameter family of extended Schwinger-Keldysh contours
We will assume that W (t) and V (0) are hermitian from here on. We formally consider the following regularization of the commutator square of Eq. (1.3): with σ ∈ [0, 1/4]. First, we note that for 0 ≤ α ≤ 1, C(t; β) (α,0) is positive definite and for α = {0, 1}, σ = 0 we recover the unregulated thermal commutator squared in the thermal state. Expanding the terms in C(t; β) (α,σ) gives Eq. (1.5) The last two are conventional Schwinger-Keldysh time-ordered correlation functions (TOCs), whereas the first two are true out-of-time-ordered correlators of the type .

(2.4)
Each out-of-time ordered correlator F (t 1 , t 2 ) (α,σ) may be seen as a correlation function in the extended Schwinger-Keldysh contour. The usual choice with α = 1/2, σ = 1/4 is shown in Fig. 1-(a); the more general F (t 1 , t 2 ) (α,σ) corresponds to a more complicated contour like the one shown in Fig. 1-(b) with different separations in imaginary time between each of the branches. It is this OTOC F (t 1 , t 2 ) (α,σ) that controls the regime of exponential growth and the Lyapunov spectrum , with A(0) a finite positive number. We will now show that the same exponential time dependence and thus the same Lyapunov exponent is obtained independent of the value of σ if α = 1/2. This follows directly from the analyticity property of the function highlighted above: [3]. Suppose for the particular value σ = 0 the function F (t 1 , t 2 ) (α,0) has the exponential behavior . Substituting this into Eq. (2.4) we get For the specific choice α = 1/2 -the one that is made in almost all previous studies -the prefactor A(iβ(1−2α mod1))| α=1/2 = A(0) is the same in both cases and equal to the one computed for the α = 1/2. Thus Although the Lyapunov exponent is not affected by the deformation parametrized by σ away from (α, σ) = ( 1 2 , 0), we do see that the prefactor of the exponential depends on the σ-deformation of the contour. Therefore, similarly to a Wightman function in Schwinger-Keldysh theory, the full commutator square C(t; β) ( 1 2 ,σ) cannot be an observable measurable in an experiment, even though it may contain physical information.
We also point out that the dependence of the prefactor on the contour seems to be in tension with the recent attempts to associate maximal chaos, defined as maximal Figure 1.
(a) Extended Schwinger-Keldysh contour corresponding to Tr ρ Lyapunov exponent λ = 2π/β, to destructive interference of the commutator square [14,15]. The destructive interference refers to the fact that, if the decoherence factor equals cos(λβ/4), it vanishes for maximal chaos λ = 2π/β. This implies that for maximal chaos the exponential time-dependence should be absent in the symmetric commutator square. Our derivation shows, however, that this is an artefact of the analytical continuation. In our case, the decoherence factor of commutator square of Eq. (2.6) is cos((1/4 − σ)λβ), which does not vanish for maximal chaos λ = 2π/β, provided 0 < σ ≤ 1/4. This casts doubts on how universal the relation between maximal chaos destructive interference may be.
Moreover, it has also been suggested that in SYK the prefactor of the OTOCs A cos(λβ/4), where A = βJ/N , is an observable which is finite at zero temperature [16]. However, as we have shown above this quantity is contour-dependent and therefore, it is not manifestly a physical observable.
To summarize, since the commutator square depends on the contour it is not clear whether the regularised commutator square is actually an observable. Another possibility may be that not all regularizations of the commutator square are physically allowed and one value of σ is preferred. For the specific deformation parametrized by σ, we could not find an argument for such case.

The α-contour
Starting from α = 1/2, the parameter σ affects only the decoherence factor of the commutator square but leaves the Lyapunov spectrum invariant. There is therefore a possibility that the Lyapunov spectrum as defined through the OTOC does measure a physical quantity. We set σ = 0 from here on and now explore its dependence on the other contour parameter α which fixes the distance between the forwards branches as shown in We have already seen that, for α = 1 2 , different choices of σ cannot be related by analytic continuation. Neither can C(t; β) (α,0) and C(t; β) (α ,0) be related to each other by analytic continuation. In other words, the distance in imaginary time between the forwards branches cannot be compensated by analytic continuation of time. This may be seen explicitly by rewriting the OTOCs in C(t; β) (α,0) as follows where we have chosen to compare to the standard contour with ρ 1/4 separation. The differences between the complexified times, t 1 + iβ(α − 1 2 ), t 2 , t 3 + iβ(α − 1 4 ) and t 4 +iβ/4 in Eq. (2.8), no longer vanish in the analytically continued OTOCs and this prevents relating one Lyapunov exponent to another. In particular, the imaginarytime separation between the two operators V (0) in both G α and H α depends on α. The standard choice, F (t) ( 1 2 , 1 4 ) , which is the building block used to derive the bound on the Lyapunov exponent [3], is computed on a contour where the separation is β/2 and α = 1/2. Therefore, G α and H α cannot be related to F (t) ( 1 2 , 1 4 ) by a simple analytic continuation whenever α = 1/2 and we have to study the behavior of these OTOCs separately.

OTOCS and physical observables in SK formalism
As one may extrapolate from the previous section, the OTOC and its Lyapunov spectrum will in general depend on the Schwinger-Keldysh contour on which it is computed. At first, this result may be surprising because, in standard Schwinger-Keldysh, it is known that physical Green's functions are independent of the contour due to energy conservation [6,7]. Indeed, since the doubling of the contour is an artificial mathematical convenience, a priori only correlation functions with external insertions on a single branch should considered physical, e.g., where we indicated with (i) the branch where each operator is inserted. With this definition, the fact that the correlation functions do not depend on the contour is a simple diagrammatic proof. We restate it here for the sake of clarity; it can be found in [6,7]. By inspecting the SK effective action, we know that the interaction vertices are of the form (2.10) Consequently, in the diagrammatic expansion each vertex is either of type 1 or of type 2. The external legs of the vertices are connected to each other or to external operator insertions with the propagators Without loss of generality, we focus on the simple lowest order 1PI diagram with n operators inserted the branch 1 and only one n-point vertex: Clearly if the vertex is of type 1, there is no contour dependence in the diagram. When the vertex is of type 2, as in Fig.3, we need to use a Wightman function. For a general contour where the forward and backward branches are separated by ρ α this is one of the Wightman functions 2 By Fourier transforming the time direction, using ρ αÔ (t)ρ −α =Ô(t+iα) and Fourier transforming back, one readily derives that  At lowest order, there is a single n-point vertex on branch 2. Contracting each of the legs of the vertex with the external operators on branch 1, and by using (2.13), this means that the relation between correlation function on different contours is Because of energy conservation at the vertex, i=1,..,n k i = 0, the overall factor vanishes and this proofs the contour independence of these types of diagrams.
However, if one of the external legs is in the branch 2, see Fig. 4, it is easy to see that now one of the Green's function no longer depends on the separation α at all, so the global factor in the n point function does not simplify anymore. The simplest example of this is the Wightman function itself. There is no vertex, but we have already shown that G α 12 = G 0 12 above in Eq. (2.12). Extending to an n-point correlation functions with a single n-point vertex, one has but now the exponent in the prefactor i=2,..,n k 0 i = −k 0 1 = 0. It is not difficult to see that the simple proof shown above extends to any diagrams. Indeed, given any diagram of the expansion, it is sufficient to divide it in subdiagrams and to use the momentum conservation in each vertex.
Turning our attention back to the OTOC, by construction each insertion occurs on one of four different branches. This indicates that the OTOC will be contour dependent, similar to two-branch correlation function in Schwinger-Keldysh theory as depicted in Fig.4. If so, this does not immediately mean that the OTOC does not measure a physical quantity (in part). For example, the (bosonic) Wightman function G βα 12 (k) = e βk 0 (1 + n(k 0 ))ρ(k) depends on the contour, but still encodes a physical quantity, namely the spectral density ρ(k). Therefore, more care is needed to understand the relation between the contour-dependent OTOC and physical properties of the system.

Contour dependence of the Lyapunov spectrum in a matrix Φ 4 theory at weak coupling
We now prove by direct computation that the OTOC indeed depends in detail on the contour chosed. In this section, we compute Lyapunov spectrum obtained from the commutator square C(t; β) (α,0) in a perturbative matrix field theory, which has been studied in detail for α = 1/2 in [8]. In the next section we will repeat this for the SYK model. The advantage of the perturbative field theory calculation is that the commutator square can be related to a kinetic equation encoding the microscopic dynamics [11]. From this, we will suggests that this microscopic insight argues that one specific contour, the one with α = 1/2 is the one that computes microscopic chaos.
We consider a 3+1 dimensional QFT with a Hermitian matrix field Φ ab whose Lagrangian is given by The commutator square of Eq. (2.7) in this matrix model is For t > 0, which we shall assume, the lowest order (disconnected) contribution is the product of two retarded Green's function arising from a contraction on the top two folds and the bottom folds separately; there is therefore no contour dependence. The non-trivial contribution at the next order, that can seed exponential growth, is the contribution with two Wightman functions connecting the two retarded Green's functions. For α = 1/2, this equals [8]: where the kernel R(p) is determined in terms of Wightman functions with operators separated by iβ/2:

4)
Note that it is only G 12 (k) and not G 21 (k), independent of the deformation α = 1/2 which appears inside the kernel. This choice is due to the identity G 12 (k) = G 21 (−k). We will consistently use G 12 only; this will not affect the final result. Defining a function f (ω, p), at the next order one of the contributions is and by rewriting C(ω) , one can set up a recursive Bethe-Salpeter equation to determine f (ω, p) and hence C(ω) ( 1 2 ,0) to all orders. Since we are interested in the late-time exponential growth, we focus on the homogeneous part of the Bethe-Salpeter equation, which in the low-frequency, late time limit equals The positive eigenvalues of the rung function R(k − p) considered as a matrix in k and p form the Lyapunov spectrum characterizing the exponential growth at late times, as we will review below. It is therefore clear that the only α-contour-deformation dependence arises from the Wightman functions in the rung function. It is then straightforward to derive the contour-dependence of the OTOC. For α = 1/2, the rung function should be modified as sketched in Fig. 5. Mathematically (3.8) Again defining this will now obey the equation: with η ≡ β(α − 1/2). Note that the change in the rung function does not depend on whether it is constructed from G 12 (k) or G 21 (−k). This can be confirmed by the fact that the commutator square should obey a KMS type symmetry α → 1 − α on the doubled time contour. This follows by redefining k → ω − k and p → ω − p. The kernel R(k − p) is even in k − p as can be readily seen from its definition Eq. (3.4). The product G R (p)G R (ω − p) changes into itself, and one obtains an equation for f (ω; ω − p 0 , −p) which identical to the original equation. In the long time low frequency limit the product of retarded Green's function is dominated by a pinching pole singularity, which amounts to the following approximation [8] where the imaginary part of the two-loop (α-independent) self energy Γ p is also determined in terms of (the α = 1/2) R(k) Eq. (3.4): (3.12) To solve the Bethe-Salpeter equation, one then makes the natural ansatz Substituting this ansatz together with eq. (3.13) into eq. (3.10), we then perform the integral over p 0 . This yiels where we have explicitly used that the rung kernel is even in the energy argument: In the time domain this is an equation of the type The solutions are the eigenvectors of M pk with an exponential growth/decay in time proportional to the eigenvalue. The positive eigenvalues of M pk are the Lyapunov spectrum. This can be found numerically; the precise method used to solve this equation may be found in Appendix A. Without computation it is already clear, however, that the result will depend on the α-deformed contour, as the defining Bethe-Salpeter equation does so. The result is presented in Fig. 6. We clearly see the dependence of the two positive Lyapunov exponents on the contour. The spectrum does become contourindependent in the high-temperature limit. This follows directly from the fact that the deformation parametrized by η = β(α − 1/2) becomes negligible for small β (compared to the mass). For intermediate and small β, the Lyapunov spectrum sensitively depends on the choice of contour. As also noted already in [8], in the extreme low temperature limit βm → ∞, the Lyapunov spectrum vanishes exponentially in βm. Even though this decreases the relative dependence on the contour, the contour dependence still persists and is given by e −(β+2|η|)m , as shown in Fig. 7.

The contour dependence regulates the IR
The contour dependence of the Lyapunov spectrum explicitly exhibited above emphasizes an important point regarding the physics behind the contour deformation. One of the arguments made for deforming contour symmetrically is that the smearing of the density matrix regulates a short distance singularity by separating the local operators in imaginary time. If this were indeed what the smearing should accomplish, then (1) at any finite value of regulator η we should expect the low-temperature limit to be universal, and (2) at any finite temperature β in units of the mass m the answer for the OTOC should diverge as one removes the regulator |η| → β 2 . The result, however, shows the opposite. The high-temperature limit is universal, indicating that this is the regime that is insensitive to the regulator, and, though we do not compute the full OTOC, the Lyapunov spectrum at fixed βm stays finite for any value of regulator. This argues strongly that the contourdeformation regulates the IR rather than the UV. This in fact agrees with Schwinger-Keldysh theory. There the "contour-deformation" is the introduction of temperature itself, and this is a well-known IR regulator.
For results in the literature in perturbative QFTs this diametrically opposite interpretation of the contour deformation has little effect. As in e.g [8,11,17,18] usually the focus is on the universal high temperature regime. However, for the SYK model, the focus has often been on the emergent regime at low temperatures. There, this realization that the contour deformation regulates the IR, may imply that the results are in fact regulator dependent and do not reflect physical information about the true dynamics. Before we turn to this in section 4, we first try to address how to obtain the physical information about the true chaos/scrambling dynamics at low temperatures.

Kinetic theory interpretation of the α-deformed OTOC
IR regulators often encode real physical circumstances. The correct question to ask therefore is which contour properly reflects physical information of microscopic chaos. In this section we will argue that this can be decided by interpreting the result of the previous section in terms of a kinetic theory for many body chaos derived in [11]. There, the authors showed that the computation of the α = 1/2, σ = 1/4 OTOC is equivalent to a Boltzmann-like equation that tracks the time evolution of the gross energy exchange.
We briefly review this result. The standard Boltzmann equation describes the time evolution of the single-particle distribution function f (t, r, p), 3 parametrizing the deviation of the single-particle distribution function from its equilibrium value:  20) where L(p, l) represents the collision integral. L(p, l) contains two contributions, namely the gain term R ∧ (p, l), counting increase of the density of the phase-space cell, and the loss term R ∨ (p, l), which accounts for scattering out of the phase-space cell. In terms of these two contribution, the Boltzmann equation is As shown in [11], the Bethe-Salpeter equation of the symmetrised commutator square C(t; β) ( 1 2 ,0) is equivalent to considering a Boltzman-like equation where the sign of the contribution of the true loss term is changed, so that we account for a gross exchange rather than a net exchange. More precisely, the gross exchange is given by where E[E p ] = 1/ sinh(E p β/2) and the extra factor Γ l , the self-energy due to the thermal environment, is present to avoid over-counting. It can be understood as follows: R ∨ T (p, l) ≡ R ∨ (p, l) − 2Γ l δ(p − l) counts the changes in the particle number f (t, p) due only to processes with p = l. Therefore, changing the sign of R ∨ in Eq. (3.21) would over-count the contribution from the bath. If one changes only the sign of the true loss term R ∨ T (p, l), the gross exchange is exactly given by R ∧ (p, l) + R ∨ (p, l) − 4Γ l δ(p − l) [11]. The eigenvalues of the integral operator (3.22) are equivalent to those measuring the exponential growth rate of the OTOC, and thus give the Lyapunov spectrum of the theory. As the α-deformation only changes the rung function in the Bethe-Salpeter equation, resulting in result (3.15), it is immediately recognized that the kinetic equation encoding the late time behavior of these families of OTOC is modified as follows The kinetic equation equivalent of the contour-dependent commutator square gives us a direct physical interpretation of what is computed, as we understand each term as loss, gain and self-energy terms in the microscopic dynamics. The explicit η = β(α − 1/2) dependence in Eq. (3.23) shows that the different contours in the α-family have a different physical origin. While for η = 0 (symmetric regularization) both the gain and loss processes are weighted equally, for other contours η = 0, their relative weight is different. For none of these values does the kinetic equation have an obvious natural physical interpretation in terms of gross, net or otherwise simple exchange dynamics.
On the other hand, the gross exchange equation has been put forward independently already a long time ago as a measure of microscopic classical chaos [13]. This conclusion from the weakly coupled field theory computation above therefore strongly suggests that, in order to probe dynamical many-body chaos in QFT, the correct choice for the out-of-time correlation function is the symmetrically regularized choice with η = 0. Fortuitously, this is the one that has predominated all the calculations in the literature, including the derivation of the MSS bound on chaos [3]. It also means that the naive thermal expectation value of the commutator square Tr[ρ[W (t), V ] 2 ] does not measure microscopic quantum chaos. One is therefore left with the reversed question: how does one justify from first principles the symmetrically regularized commutator square as a measure of quantum chaos. We will return to this question in the last section. First, we will consider the same question of contour-dependence of the commutator square, its Lyapunov spectrum and an argument in favor of symmetric regularization in a strongly coupled field theory -the SYK model.

Contour dependence of the Lyapunov exponent in the SYK model
One of the research directions where the commutator square has had important impact is in the emergent strongly coupled low energy regime of the Sachdev-Ye-Kitaev model. The exponential growth of the symmetrically regularized commutator square saturates the MSS bound on chaos λ L ≤ 2πT and this has given great impetus to the notion that the SYK model provides a microscopic theory for AdS black holes. Now that we know that the commutator square and its Lyapunov spectrum depend on the way the contour is regulated, the natural question arises how this affects the insights in the SYK model. We shall first show that indeed also in the SYK model the Lyapunov spectrum is contour regularization dependent. As before there is a universal contour independent high temperature result, but the contour dependence is strong in the low temperature IR -the very regime which of most interest in the literature. In fact, for any contour other than the symmetric one, the Lyapunov exponent violates the MSS bound. Granted, this is not a fair comparison as the MSS bound is specifically computed for the symmetrically regularized contour. Nevertheless, it re-emphasizes the question, which regularization computes the physical microscopic scrambling in the system. In weakly coupled QFTs the kinetic equation interpretation suggested that fortuitously it is the symmetric configuration that does so. It turns that also for the SYK model there is a different computational route to extract the true scrambling rate: and again this happens to agree again with the symmetric regularization.
The SYK Hamiltonian with q/2-body interactions is where χ i are Majorana fermions so {χ i , χ j } = δ ij and the coupling J i 1 ,i 2 ,...,iq is a Gaussian-distributed random variable with zero average and diagonal (i.e. for each [19]. The fermionic two-point function G(τ ) = − T χ(τ )χ(0) satisfies the following averaged Dyson equation in the large-N limit [19]: with ω n = (2π/β)(n + 1/2), G n ≡ G(iω n ) and Σ n ≡ Σ(iω n ). In the same way as for weakly coupled QFT the symmetrical contour regularized commutator square C(t; β) ( 1 2 ,0) satisfies a Bethe-Salpeter equation. In the large-N limit, for arbitrary coupling, this BS-equation is [19]: where G R and G W are the retarded and Wightman two-point functions obtained by analytic continuing the Dyson equation Eq. (4.2) to real time and solving these equations numerically with an iterative procedure [20]. 4 In slight contrast to the perturbative QFT approach, in the SYK literature one solves Eq.(4.3) by making the explicit ansatz F (t 1 , t 2 ) = e λ L (t 1 +t 2 )/2 f (t 12 ) and rewriting it as an integral eigenvalue equation in frequency space: One then (numerically or analytically) searches for which value of λ L the kernel has an eigenvector with eigenvalue 1 [19].
As in the matrix model of the Sec. 3, the only place the contour regularization shows up is in the Wightman functions. 5 Instead of parametrizing with respect to the α = 0 Wightman function, let us parametrize with respect to the α = 1/2 Wightman function: (4.5) The Bethe-Salpeter equation (4.3) for the commutator square in frequency space for arbitrary α-deformed contour is then the same as before, but with a modified kernel g lr (ω):f We evaluate the modification in the kernelg lr (ω) compared to the original kernel g lr (ω) by using the convolution of the Wightman functions: and substituting G η (ω) = e ωη G W (ω) in each term inside the integral: Therefore, Eq. (4.6) reduces tõ Up to the inclusion of the exponential ansatz, which makes this equation dependent on λ, this is essentially the same equation as the perturbative QFT Eq. (3.10). Numerically determining the Lyapunov exponent for which the kernel has eigenvalue 1, we show in Fig. 8 the result as a function of the coupling for different choices of the regularization α of the OTOC G α (t) (2.8). As expected, one sees also for the SYK model in the large N -limit a clear dependence of the Lyapunov spectrum on the contour parameter α. Two aspects are noteworthy. Also as predicted, the regularization dependence is strongest in the IR, but this is precisely the regime of most interest in SYK studies. Furthermore, the dependence is strong enough to go above the bound λ ≤ 2π/β proved for the symmetric regularization η = 0 [3], an observation also made by [5] for a system of 2D fermions with quenched disorder and η = β/2 (in our notation). Refocussing on the question which of the regularizations then encodes microscopic process of scrambling/chaos, we now consider the two regimes of the SYK model where some analytical control is possible: the strong coupling limit βJ 1 and the large-q limit.

Study of the OTOC in SYK in the strongly coupled limit: conformal limit analysis
In the strongly-coupled regime βJ 1 of the SYK model, where conformal symmetry emerges asymptotically, the OTOC may also be computed analytically by studying the spectrum of the Casimir operator. More specifcally, for βJ 1 the eigenvectors of the Casimir operator, with eigenvalue h(h − 1), are also eigenvectors of the Euclidean kernel of the Bethe-Salpeter equation [19]. In this regime, the kernel of the Bethe-Salpeter equation is: where the eigenvalues of K c depend on q and h. Moreover, the allowed values of h are constrained, because the Bethe-Salpeter equation for the OTOC selects the eigenvalue unity of K c . For q = 4, the leading contribution to the OTOC turns out to be h = 2 and is given by [19]: where θ is the rescaled Euclidean time θ = τ /β. This equation must now be analytically continued to real time by choosing the operator insertions. We consider Figure 9. Extended Schwinger-Keldysh contour corresponding to the two-parameter OTOC Tr the contour shown in Fig. 9, which allows us to consider both the σ-and α-families simultaneously.
More specifically, we choose In terms of x, x , y, y , we have: In order compute Eq. (4.9) explicitly, we set x = x , sum over n and then substitute Eq. (4.10) to get: which for large t behaves as: We first note that F(t) is symmetric over α → 1 − α, as expected. Second, the long-time regime is controlled by a growth rate given by 2π/β, independent of the distance between the forward branches α. This result seemingly contradicts the numerical results shown in Fig. 8. However, this may be understood from Eq.  [19]. Therefore, the η = 0 kernel g lr (ω) is Using the identity |Γ(a + ib)| 2 = Γ(a) 2 ∞ k=0 1 1+b 2 /(a+k) 2 , one immediately sees that this kernel is strongly peaked around the origin β(ω − ω) → 0. On the other hand, changing the regularization changes the kernel by an overall factor e (α−1/2)β(ω −ω) . Thus, as the integral in Eq. (4.8) is dominated by ω ∼ ω, the dependence on the contour proportional to e η(ω −ω) ∼ 1 essentially drops out. This explains the contour independence of the strongly coupled conformal limit analysis of the SYK commutator square, and the "contour-independent" result is equivalent to the η = 0 result computed in the conventional Bethe-Salpeter equation. This supports our kinetic equation finding that the symmetric contour encodes the true physics.
From this, we conclude that the different lines in Fig. 8 should eventually converge at λβ/2π = 1, for βJ → ∞. This regime, however, is not accessible in our numerics.

Study of the OTOC in SKY in the limit of large interaction order
In the SYK model, analytical control is also possible when one increases the order of the interaction in the Hamiltonian, which is controlled by q [19]. Here we consider the calculation of the Lyapunov exponent in SYK in the large-q limit, and show that it is also contour independent and reproduces all qualitative features of the symmetric regularization of the OTOC at finite q: it is monotonic and does not violate the MSS bound, and becomes identical in the limit βJ 1. We start with the two-point function in Euclidean signature in the large-q expansion, where q controls the order of the interaction of the Hamiltonian Eq. (4.1) [19]: where g(τ ) is obtained by inserting the above ansatz in the saddle point equation for the two-point function. This gives the equation where θ = τ /β ∈ [0, 1) and J 2 = q2 1−q J 2 , and with boundary conditions g(0) = g(1) = 0. The solution of Eq. (4.14) where ν ∈ [0, 1] parametrises the flow from weak βJ ∼ 0 coupling for ν ∼ 0, to strong coupling βJ 1 for ν ∼ 1. The analytic continuation to real time: for α = 1/2, G (α) (t) gives the Wightman function with operators separated by iβ/2. Instead of working in frequency space with Eq. (4.8), we work in the time domain and use the following simplification for large q: .
Here we are able to find numerical evidence that supports the finding that the large q-approximation commutator square is independent of the contour. In Fig.  10 we show that, indeed, for q = 10, the dependence of the Lyapunov exponent on the contour parameter is weaker than for q = 4. This suggests that indeed, in the q → ∞ limit, all the regularizations give the same Lyapunov exponent equal to the one computed for the symmetric regularization with η = 0.

The Lyapunov spectrum and the Loschmidt echo
In previous sections we have seen in two different theories the inevitable regularization dependence of the commutator square. This shows that without more detailed specification one cannot directly relate this quantity to an observable that can be measured in experiments. We have also shown that the regularization dependence is dominant in the IR rather than the UV. This is analogous to Schwinger-Keldysh theory where contour dependence is related to the temperature, and the latter is a well known IR regulator. IR regularization issues are usually not solved by counterterms and renormalization. Instead they often encode physics on their own. This suggests that a way to resolve the regulator dependence is to define which member of family of "regularized" correlation functions computes a proper physical observable. Both the weakly coupled QFT result through the mapping of the commutator square to a kinetic equation and the strongly coupled SYK model through direct alternate computation of scrambling in the IR indicate that the symmetrically regularized commutator square is the correct one.
Fortuitously this the one almost exclusively studied in the literature and the one for which the MSS bound on chaos is derived. Nevertheless one would like to understand from first principles why the symmetrized contour is an appropriate physical observable. The first attempt construction in the introduction points to the thermally averaged commutator square instead. In this section we show that the symmetrized commutator square follows directly from an alternative measure of chaos, which is related to standard measurements of information spreading: the Loschmidt echo. This quantity contains not only the commutator square but also higher-order out-of time correlation functions. The Loschmidt echo and related quantities have been used in the context of quantum chaos for a long time [21][22][23][24][25]. Therefore, it is not surprising that the OTOCs may be extracted from echo spectroscopy as proposed in [26][27][28][29] and measured experimentally in [10].

Loschmidt echo
The Loschmidt echo is based on a old thought experiment trying to disprove the irreversibility inherent in Boltzmann's equations by imaging a dynamical system where at time t after t 0 = 0 one reverses all velocities and compare the resulting state at time 2t with the original state. Microscopically the answer is of course identical, but supposing one makes a tiny "erroneous" perturbation at the time when one reverses all velocities, one immediately sees that in a chaotic non-integrable system the resulting state will be exponentially different from the original state.
This thought experiment can be directly mapped to a quantum quench experiment. One evolves a quantum state forward in time for a time t, perturbs it with an instantaneous quench e iδW , evolves backward for the same time t and projects onto the original state.
For a generic initial state, the echo will have a universal late time exponential fall off independent of the type of quench W that encodes the lack of overlap between the initial and final state.
The Lyapunov exponent λ is then a property of the system characterized by its Hamiltonian H alone. The Loschmidt echo is the expectation value of a complex operator. To avoid phases one often takes the absolute value squared, which is known as the fidelity [30] The intermediate step is a well-known result from Jalabert and Pastawski [22]. A second practical step with an eye on experiment is to consider the fidelity for an ensemble of states, rather than a single state. Choosing the thermal ensemble one has Defining e −iHt ρe iHt ≡ X and e −iδW e −iHt ρe iHt e iδW = Y , the fidelity F (t) above is a specific case of the more general operator fidelity F = Tr X † Y applied to density matrices as operators. 6 Three remarks are in order. (1) In the intermediate step we used that the leading Lyapunov decay rate in t is the same when computed via Trρe −iδW (t) 2 or Trρe −iδW (t) ρe iδW (t) .
(2) Naively, as the late time Lyapunov exponent of interest is a property of the system and not of the initial state, the averaging should not matter. However, it is well known from classical dynamical systems that the late time behavior of an ensemble of classical trajectories is governed by Policott-Ruelle decay, rather than the microscopic exponential growth. Even though these are qualitatively related in weakly coupled theories, they are not quantitatively the same [11]. (3) Note both the symmetrized appearance of the density matrix, and the fact that the cumulative power of the density matrix is 2. Computed through a path-integral this implies that the periodicity in Euclidean time is twice the inverse temperature β = 2/T phys . 6 The operator fidelity is a weaker version of state fidelity encoding the notion of how close a state is to a maximally entangled one [31] or, if referring to teleportation, it quantifies the quality of the teleportation that can be achieved with the given state [32]. The state fidelity between two quantum states given by the density matrices ρ 0 and ρ 1 equals [33,34]: To connect with the commutator square, we expand to second order in δ F (t) = Tr ρ 2 + Tr ρ(−δ 2 W (t) 2 )ρ + Tr ρ(δW (t))ρδW (t) = Tr ρ 2 + δ 2 2 Tr [ρ, W (t)][ρ, W (t)] + . . . (5.6) with the difference that the density matrix itself takes the role of the operator V (0). The second time dependent term, the density-matrix commutator square, is a variant of the Wigner-Yanase-Dyson skew information.
for the symmetric value α = 1/2 [35]. Writing out the symmetric case for hermitian A, and replacing the thermal density matrix ρ with a pure state density matrix one can recognize that the WYD skew information is an extension of the variance for pure states to mixed states. If, by the same argument as above, one may assume that it is dominated by some largest eigenvalue Tr ρAρA ∼ (Tr ρA) 2 , it computes something akin to the (largest eigenvalue) variance for the operator O = ρA. In that sense it is again natural that the density matrix appears with cumulative power 2. Put differently, in computing the WYD skew information the periodicity in Euclidean time is twice the inverse temperature β = 2/T phys . However, this is not yet the commutator square we are interested in. A guess might be the case where the thermal density matrix is rotated by a small similarity transformation ρ = e iV ρ 0 e −iV . This is equivalent to an instantaneous quench by V at time t = 1. Then in the limit of small δ the late time fidelity equals F (t) = Tr ρe iHt e −iδW e −iHt ρe iHt e iδW e −iHt (5.10) The first and the last term can never give an OTOC, however, and, ignoring those, one has in the limit of small V The two terms of order δ 2 in the first line are also TOC. The terms on the second line contain the symmetric commutator square and a second term which also an OTOC but on a different contour. 7 As we know by now, generically the Lyapunov behavior of this other OTOC will be different. This is not yet the answer. Tracing the origin of Eq.(5.11), it is easy to see how the fidelity and the symmetrized commutator square are related. Eq.(5.11) follows from taking the long time limit and then taking V and δW small in the fundamental definition of the ensemble averaged fidelity -the first line of Eq. (5.4). If, however, we take the limit of V and δW small, with ρ = e iV ρ 0 e −iV , the ensemble averaged fidelity equals We now use the late time approximation, where we assume that ρ 0 [V, W (t)] is dominated by an eigenvalue Eig(ρ 0 [V, W (t)]) ∼ e 1 2 (λ+iφ)t . In that limit the middle two terms give a strongly oscillatory contribution, which is hard to measure. We therefore ignore it. As to the last term in Eq.(5.12), there the late time limit allows us to make again the approximation We recognize precisely the symmetrized commutator square with one already noted difference. The cumulative power of the density matrix is 2. This implies that the connection between the periodicity in Euclidean time and the physical temperature differs with a factor two compared to what the naive smearing procedured assumes: β = 2/T phys . In particular this means the proper MSS bound on chaos should read λ ≤ πT phys . The above is a strong argument that the natural observable which measures the symmetrized commutator square is the Loschmidt echo in the limit of small quenches first and late time subsequent with the sublety that β = 2/T phys .

Conclusion
In this article we have explored the role of the regularization scheme of the commutator square and of the OTOC. Quantum chaotic systems may display an exponential growth parametrized by a quantum Lyapunov exponent which is bounded by above λ ≤ 2πk B T / [3]. The proof of this bound involves regularising the OTOC by thermally spreading the operators. Purportedly, this is done to regulate short distance singularities and any physical property of a system should be independent of the short distance regularization scheme.
Here, we have shown that for those regularizations consisting on a contour with a iβ/2 separation between the forward branches, shifting the backwards branches induces a change in the decoherence factor, defined as the prefactor of the sum of the OTOCs [14]. Therefore, the decoherence factor cannot be a physical quantity as previously suggested. On the other hand, the Lyapunov exponent is the same for all of these contours, suggesting that indeed it may be measurable. However, we have then shown that for a different choice of contours, where the separation between the forward branches is changed, the Lyapunov spectrum also depends on the contour chosen. While the contour dependence of the commutator square has been mostly overlooked in the literature, it is not surprising that this is the case. Similarly to the Wightman function, the commutator square involves operators inserted on forward and backward branches of the Schwinger-Keldysh contour, and so there is no reason to expect that it should be a physical quantity. Therefore, it is important to know how to extract physical information from it, in the same way that the spectral density, a physical quantity, may be obtained from the Wightman function, even though the Wightman function itself is not physical. The one notable exception in the literture is [5]. These authors studied many body chaos in a weakly interacting 2D system of fermions with quenched disorder and computed the Lyapunov exponent both for the unregularized η = 1/2 case and the symmetrically regularized one, η = 0 (in our notation). They indeed found that the two results disagree, and pointed out the regulator dependence of the OTOC. The conclusion that they drew is that, in the model considered, the only special feature of the symmetrically regularized OTOC is a particular cancellations of divergencies in the computation, but the physical meaning behind this correlator remained obscure.
We have done a more thorough analysis showing the regulator dependence of the OTOC for two paradigmatic models, a weakly coupled φ 4 matrix boson (at any N ) and the strongly coupled SYK model. By comparing to ordinary Schwinger-Keldysh theory, we provide a simple diagrammatic proof regarding the reason why the choice of the contour affects the OTOC. Moreover, we exhibited how for the SYK model, which has been extensively studied over the last years, the results near the low energy conformal regime of most interest strongly depend on the contour to the extent that the other regularizations violate the bound.
These detailed studies allow us to recognize that the regulator dependence is an IR issue, and not an alleviation of purported UV singularities. This means one has to take more care in understanding the role of the regulator as it may contain physical information. One crucial insight of our paper is to recognise the special physical meaning of the symmetrically regularised OTOC, by means of kinetic theory [11]. That the fact that the bound on chaos holds for the physically meaningful definition of OTOC is remarkable and open new directions on possible still unknown dynamical constraints that the bound can impose. We were able to confirm the special nature of the symmetrically regularized contour even in the strongly coupled limit of the SYK by using two other computational routes instead of the standard Bethe-Salpeter resummation.
This does then raise the question which simple observable naturally gives rise to such a symmetric insertion of a thermal density matrix. We proposed a simple observable, related to the operator fidelity, which contains information beyond the commutator square and can be measured experimentally using echo spectroscopy. The corollary of using this observable to define the OTOC is that is based on a double insertion of density matrices, i.e. the periodicity in Euclidean time is twice the inverse temperature. From this point of view the bound on chaos should read λ ≤ πk B T phys / .
Overall, our results pose the question on the usefulness of the commutator square to probe quantum chaos. The contour dependence of the commutator square and of the Lyapunov spectrum extracted from it, cast doubts on whether the commutator square is physical and how physical information should be extracted from it. However, even though a natural way to define chaotic quantum system is that in which the OTOC displays an exponential growth, this growth regime actually clashes with the other notion of a quantum chaotic theory that it should display random matrix behaviour. In the SYK model, even though one has exponential growth at shorter times similar to classical weakly interacting chaos, spectral properties such as the spectral form factor, are similar to that of random matrix theory for times of order of N log(N ) and larger [37]. This suggests that the model becomes truly quantum chaotic after this time-scale. A gorgeous example of true quantum chaos embodied by random matrix behaviour has been observed on the kicked Ising spin-1/2 chain for much shorter timescales [38]. There is no exponential growth in the OTOC in this model, which challenges the notion of how quantum chaos and especially maximal chaos should be defined.

A Numerical calculation in matrix model
In this appendix we outline the simplifications used to solve numerically the Bethe-Salpeter equation Eq. (3.15). Following [8], we define P = |p|, K = |k|, y = |k − p| (A.1) and use symmetry to simplify the momentum integral: Rewriting Eq. (3.15) in the time domain and replacing the momentum integral we arrive at the simplified version of the Bethe-Salpeter equation, which we solve numerically following the strategy described in [8]: where I + (P, K) ≡ K (2π) 2 P 4E P E K K+P |K−P | dyyR(E + , y) = 3λ 2 K (2π) 3 P E P E K sinh(βE + /2) K+P |K−P | dy log sinh x +