Power-law out of time order correlation functions in the SYK model

We evaluate the finite temperature partition sum and correlation functions of the Sachdev-Ye-Kitaev (SYK) model. Starting from a recently proposed mapping of the SYK model onto Liouville quantum mechanics, we obtain our results by exact integration over conformal Goldstone modes reparameterizing physical time. Perhaps, the least expected result of our analysis is that at time scales proportional to the number of particles the out of time order correlation function crosses over from a regime of exponential decay to a universal $t^{-6}$ power-law behavior.


Introduction
The Sachdev-Ye-Kitaev (SYK) model [1,2] is a system of N (Majorana) fermions, χ i with i = 1, . . . , N, subject to a four-fermion interaction with Gaussian distributed random matrix elements J i jkl with zero mean and a variance given by |J i jkl | 2 = 6J 2 /N 3 . The seeming simplicity of this model is deceptive. At low excitation energies it exhibits an asymptotically exact conformal symmetry [2,3,4,5] which manifests itself in the infinite-dimensional freedom to re-parameterize time, t → f (t), in the description of long time correlations. This 'nearly conformal symmetry' (NCFT) [4,6] makes the model a candidate holographic shadow of some two-dimensional bulk. The potential realization of a holographic principle of lowest possible dimension has triggered a surge of research activity on the SYK model and its quantum dynamics [7,8,9,10,11,12,13,14,15,16,17,18,19,20].
At the same time, the existence of a large N parameter within an 'infinite range' interaction framework make the system amenable to mean-field approaches. It turns out that at the meanfield level the infinite dimensional conformal symmetry gets broken by the interaction self-energy down to the conformal group SL(2, R) of rational transformations, t → at+b ct+d , ad − bc = 1. This leads to a classic symmetry breaking scenario and the emergence of Goldstone modes whose fluctuations become unhampered in the long time limit where the explicit symmetry breaking (represented by the time derivative ∂ t present in the system's action) becomes negligible. The situation bears similarity to that in a magnet, with the important difference that the dimension of the Goldstone mode manifold is infinite, while the spatial dimension is zero. This means that Goldstone mode fluctuations are enhanced by the Mermin-Wagner principle and expected to become virulent in the long time limit where the explicit symmetry breaking vanishes.
In Ref. [21] we introduced a method to isolate these fluctuations and perform a full integration over the Goldstone mode manifold. The idea is to introduce an exponential reparameterization f (t) = exp(φ(t)) whereupon the Goldstone mode integral assumes the form of a path integral over φ, with an (approximately) time local action coinciding with that of so-called Liouville quantum mechanics [22]. Standard methods of quantum mechanics are then applicable to perform the integration, including the regime where the φ-fluctuations become large (but those of a canonically conjugate 'momentum', k, are small in exchange). We demonstrated emergence of a new time scale At long times t M, Goldstone mode fluctuations are strong and act to restore the broken symmetry (much like in a low dimensional magnet at weak external field rotational symmetry is restored by unbounded fluctuations in the magnetization.) For example, while the mean-field Green function [1,2,4], G( ) ∼ | | −1/2 has a divergent amplitude at low energies, the inclusion of Goldstone mode fluctuations shows [21] that G( ) ∼ | | 1/2 for energies | | M −1 , i.e. symmetry restoration supresses the mean-field propagator (the latter playing the role of an average magnetization in the magnetic metaphor).
In Ref. [21] we analyzed the effect of Goldstone modes within a zero temperature framework. However, the majority of observables describing the dynamical behavior of the SYK-modelfour point functions in general, and out-of-time-order (OTO) correlation functions in particular [6,23,24] variants of the partition sum [3,24,25], etc. -are formulated as quantum thermal averages and require a finite temperature formalism. That this generalization is not entirely innocent is indicated by the observation that in the exponential degradation [24,6], ∼ exp(2πT t), of OTO correlations at small times, temperature, T , itself features as the relevant rate.
Below we show that finite temperatures affect the effective quantum mechanics describing the Goldstone mode fluctuations via the appearance of an exponential potential adding to the native potential of Liouvilian quantum mechanics. This contribution strongly affects the Goldstone mode integral. To be specific, we consider the OTO correlation function where Z = tr e −βĤ is the partition sum. This expression, which differs from the standard definition of a thermal correlation function by the symmetrization of thermal weights [24,6], has become an established tool for the diagnostics of correlations in the model. While previous work identified regimes of exponential decay of OTO correlations functions, here we show that Goldstone mode fluctuations are responsible for the formation of power laws at large times and/or low temperatures. Our main findings are summarized as follows. For high temperatures, T M −1 , we need to discriminate between short and large times, t M and t M, respectively. In the short time regime, Goldstone mode fluctuations are weak, and one identifies [6,24] two different regimes characterized by the exponential loss of correlations, 2 Figure 1: Results for the OTO correlation function. Top: At high temperatures, T > M −1 and large times, t > 2πM, the function crosses over from exponential to power-law decay with an exponent t −6 . Bottom: at low temperatures, T < M −1 the function is nowhere exponential. At large times t > T −1 > M −1 it again shows t −6 power-law behavior. The inset shows the parametric extension of the four regimes in a t − T plane.
depending on whether t < t E , or t > t E , where plays the role of an Ehrenfest time [23] in the problem. However, at large times, t M, fluctuations are strong and generate universal power law scaling, F(t) ∼ t −6 . For small temperatures, T M −1 , Goldstone mode fluctuations affect the picture throughout the entire domain, no exponential regimes are found, and the power law, F(t) ∼ t −6 , holds for all times t T −1 . The quantitative summary of these statements reads as (cf. Fig. 1 where the four distinctive regimes are indicated as 1, 2, 3, 4, respectively.) where the algebraic profiles extend previously obtained results with exponential behavior [4] to the regime of long times/low temperatures. These correlations appear as a robust consequence of 3 the Liouville quantum mechanics which effectively governs the long time behavior of the system. Though we used 4-fermion model, Eq. (1), as an example -the qualitative features, including t −6 law, hold for an arbitrary number q ≥ 4 of fermions in the interaction term. This universality is based on the 3/2 power-law decay of two-time-points functions in the Liouville model [26,21]. The derivation of the above correlation laws also implies independent validation of various results that have been obtained before on more phenomenological grounds (including by reference to principles of holography), or by different analytic methods. Notably, we find that the partition sum at small temperatures T J scales as in agreement with the results of Ref. [10]. This formula is directly related to the many-body density of states (DoS), which we find for energies E J behaves as where the energy is counted from the ground state. This result was obtained previously in the limit of large number of interacting fermions, q ∼ √ N, [10] and by statistical analysis of moments of the random interaction operator [13]. In our present analysis, the DoS reflects the energy stored in Goldstone mode fluctuations. We also note that at high temperatures, T 1/M, or energies, E 1/M, the constant M in Eqs. (6) and (7) acquires a week logarithmic dependence on either T or E, respectively (see more comments on this in section 3).
Interestingly, the result (6) for the partition sum entails an heuristic explanation for the powerlaw formation in the OTO correlation function (5) at long times/low temperatures, t, T −1 < M. To see this, let us insert four spectral resolutions of (realization specific) many body eigenstates |m in (3) to obtain the 'Lehmann representation' Now, let us assume that the Majorana matrix elements m|χ i, j |n and the eigenenergies, n , are statistically independent. This assumption is of course grossly ad-hoc. However, statistical wave function/eigenvalue independence is a hallmark of random matrix models, and in view of the evidence for ergodic chaotic behavior shown by the SYK model at low energies [20,13] may contain some truth. If so, and if at the characteristic energy differences | m − m | ∼ t −1 ∼ M −1 much larger than the many-body level spacing ∼ e −N/2 statistical correlations between the levels of different spectral sums can be neglected, the correlation function simplifies to where Z(s) = n e −s n = dE ρ(E) e −sE with Re(s) > 0 is the Laplace transform of the average many-body DoS, Eq. (7), which at s = β is also identical to the physical partition sum. The large time behavior of such Laplace transform is dominated by branching points of the DoS function. Since at energies E < M −1 DoS exhibits non-analytic behavior as ρ(E) ∼ E 1/2 , it leads to Z(s) ∼ s −3/2 at |s| M. This is exactly the 3/2 law, observed in the long-time behavior of two-point functions [21]. It in turn implies F(t) ∼ t −6 as announced in Eq. (5). 4 Of course, the construction above contains several ad-hoc assumptions and isn't trustworthy in its own right. The rest of the paper is devoted to a first principle validation of Eq. (5) from the low energy theory of SYK Goldstone mode fluctuations. We will start in section 2 where we review the path integral description of the model and the emergence of conformal soft modes. In section 3 we construct the finite temperature Liouvillian soft mode action. In section 4 the short time fluctuation behavior of this effective theory will be studied by stationary phase methods on the example of the partition sum. Continuing with this quantity the analytical machinery required for studying strong Goldstone mode fluctuations are introduced in section 5. The following four sections then define the core of the paper in which the OTO correlation function is addressed. We start by setting up a path integral representation of this quantity in section 6. Its short, intermediate, and long time behavior are then addressed in sections 7, 8, and 9, respectively, before we conclude in section 10. Several Appendices provide details on technical calculations.

Preliminaries
Our starting point is the effective action describing the system [3,4] after the averaging over disorder and integration over Grassmann-fermion degrees of freedom have been performed (see Appendix A) Here, G ab τ,τ and Σ ab τ,τ are matrix fields carrying an index structure, a, b = 1, . . . , n in a space of n replicas. The first term of the action resembles the free energy of a fermion system and identifies Σ as a dynamical self energy. The second term reflects that the average over disorder, Ĥ 2 is of eighth order in Majorana operators. This is amounts to four Green functions, which here play a role of dynamical matrix fields, too. The third term expresses the conjugacy of the two fields Σ, G.
The action (10) exhibits an exceptionally high level of symmetry, provided ∂ τ term is neglected: it is then invariant under reparameterizations of time, τ → f (τ), where f (τ) may be any invertible and differentiable function. (Invertibility requires globally monotonicity and without loss of generality we assume f (τ) > 0 throughout. The condition of differentiability may be sacrificed at isolated points if necessary.) This means that S possesses an infinite dimensional symmetry group whose generators are the Virasoro generators of infinitesimal reparameterization transformations, f .
The interpretations of the fields G and Σ become more tangible at the mean-field level, where the stationary phase equations [3] assume the form of the self-consistent Dyson equation The first equation here is a matrix equation while in the second the cube operation acts on each matrix element separately, i.e. [G] 3 ≡ G ab τ,τ 3 . These equations can be solved in the long time limit τ − τ 1/J, where ∂ τ may be neglected. A configuration, G ∞ , Σ ∞ solving the equations at zero temperature, β → ∞, reads as where b = (4π) −1/4 . However, this solution is not unique [2,5,4]. Under the action of the symmetry operation it transforms as i.e. we encounter the classic scenario where a symmetry is broken at the mean field level, and a Goldstone mode manifold emerges. One may verify that the sole transformations leaving the saddle points invariant are the conformal reparameterizations, τ → aτ+b cτ+d , ad − bc = 1, i.e. the Goldstone mode manifold is identified as the group of time-reparameterizations modulo the unbroken SL(2, R)-transformations.
The transformation (13) provides a way to generalize the saddle-point solution from one defined on the entire real axis to a finite-temperature solution with support on [−β/2, β/2] [27]: the reparameterization τ → g(τ) := tan πτ/β , provides an invertible map from the open imaginary time interval [−β/2, β/2] onto the entire axis. According to Eq. (13) the corresponding Green function reads and similarly for Σ β (τ 1 − τ 2 ). This defines a periodic and time translationally invariant solution of the finite temperature saddle-point equation (11). The group property of the symmetry implies that reparameterization of compactified time τ → f (τ), mapping the interval [−β/2, β/2] onto itself and obeying the periodicity constraint generate a family of finite temperature solutions G β ([ f ], τ, τ ), periodic on the time interval but lacking translational invariance (i.e. the functions G β [ f ] depend two time arguments separately and not just the difference).

Soft mode integration
The low energy properties of the model are described by a functional integral over the softmode manifold parameterized by the functions f identified in the previous section. In Ref. [21] we showed that at zero temperature the emerging integrals can be understood as path integrals of Liouville quantum mechanics. In the following, we discuss the non-trivial generalization to finite temperatures to then apply it to the computation of observables.
Rather than integrating over functions f : The function h then necessarily has a singularity at some τ * ∈ [−β/2, β/2], where f (τ * ) = ±β/2 (cf. Fig. 2.) A shift of the time axis, which is part of the SL(2, R) invariance manifold, can be applied to shift the singularity to the boundaries of the interval τ * = ±β/2.
in a more explicit representation. Notice that, although we effectively returned to the zero temperature Green function, the time arguments lie in the finite temperature interval τ, , formulated in terms of reparameterizations, f , of the finite temperature soft mode manifold may be transformed into the zero tempera- Expectation values of such observables are obtained by integration over reparameterizations, where we omitted the subscript ∞ for brevity. Here, Requiring invariance under the SL(2, R)-reparameterization group and assuming locality in time, one may argue [4] that this action must be the integrated Schwarzian derivative where 2 and M is a coupling constant of dimensionality [time] 1 . On dimensional grounds, this constant determines the time scale, t * ∼ M, above which the reparameterization fluctuations become strong. The value of this constant can be fixed by a more explicit construction [21] which first translates the starting action (10) into a non-local f -action and then notes that the collapse to an effectively local action takes place below a certain energy scale which has to be determined by a self consistency-analysis [21].
In more concrete terms, one defines M as the running coupling constant of the Schwarzian theory, M(E), which has a weak logarithmic dependence on a typical energy involved when the latter is sufficiently high, E > M(E), namely This logarithmic renormalization stops at the scale ∆ defined self-consistently via ∆ = 1/M(∆), which with log accuracy is resolved as the Goldstone mode fluctuations proliferate and M(E) saturates to the constant value M ≡ 1/∆ as stated in equation (2). Note, that a logarithmic dependence of M(E) is the price to pay to reduce the intrinsically non-local action for reparametrizations to the local Schwarzian form. In practical calculations one should take M at the scale E = min(T, 1/t), where t is a time variable in the correlation function. Before turning to the actual integration procedure in Eq. (18), we impose one more change of variables and introduce the degree of freedom, φ(τ), of Liouville quantum mechanics as where the positivity h (τ) > 0 guarantees that φ(τ) is real. The key advantage of this parameterization is the flatness of its integration measure, µ[h] Dh = Dφ [21]. However, attention must be payed to the singularity of all admissible h-functions at τ = ±β/2. It translates into a diverging integral We find it convenient to regularize this divergence by introducing a finite but long time t H such that and take a limit t H → ∞ at the end of the calculation. Trajectories φ(τ) satisfying this condition thus become regularized at the boundaries of the time interval, φ(±β/2) = φ 0 with φ 0 → +∞ as t H → ∞.
In the φ-language the soft mode action reads as and expectation values of observables O β assume the form where γ is a Lagrange multiplier enforcing the condition (22), and Z β = Dφ exp(−S [φ]) is the partition sum. Notice that the full derivative φ cannot be ignored since, due to singularity discussed above, φ (±β/2) → ±∞. It is this term which neutralizes the formally infinite constant γt H . The same singularity also shows in the boundary conditions φ(±β/2) = φ 0 . Our treatment below establishes a balance between the parameters t H , φ 0 , and γ such that the final answers do not depend on regularization specifics.

Saddle point calculation
At sufficiently high temperature, and/or sufficiently short time intervals τ − τ the functional integral (24) may be evaluated in a saddle point approximation. This saddle point approach is stabilized by the largeness of M compared to all other time scales in the problem. To expose this parameter in the best possible way we introduce a rescaled variable as where we are simplifying the notation by temporarily setting β ≡ π. In the final results the βdependence can then be re-introduced by scaling all quantities with dimension of time as τ → τπ/β. Expressed in the new language, the action assumes the form the constraint becomes and the equation of motion reads as ϕ (τ) = 2e ϕ(τ) . In the limit t H → +∞ this is solved bȳ The meaning of this solution becomes apparent when we formulate it in the language of the h-representation, , where g is the tan-function of Eq. (14). It is straightforward to check reduces to the familiar [27] translationally invariant finite temperature solution of the Dyson equation (15). Finite values of t H regularize the infinity in the boundary conditions ofφ ∞ . This regularization can be achieved by shifting the argument of the cos in the solution as, For small λ we find that the integral in (27) evaluates to 4/(πλ), and this fixes the shift parameter as For this configuration, we find a boundary value, exp(φ(±π/2)) ≡ exp(ϕ 0 ) = (2/πλ) 2 . Now is a good time to discuss the status of the Lagrange multiplier, γ. This parameter was introduced to enforce the constraint, Eq. (27). In conventional variational calculus, the value of a Lagrange multiplier is fixed by requiring that the constraint is satisfied on the variational configuration (φ). Presently, however, we have one more free variable in the picture, viz. the diverging boundary values ϕ 0 . We may use this freedom to make sure that for any given value of γ the constraint is satisfied (we saw above that this is achieved by the value exp(ϕ 0 ) = (2/πλ) 2 = (γt H /4M) 2 .) This gives us the freedom to fix γ at will, and the best choice is γ = J, i.e. the Finally, a straightforward calculation shows that the stationary phase action is given by independent of the regularization. Gaussian fluctuations aroundφ lead to a fluctuation determinant whose calculation is detailed in Appendix C. Multiplying this factor with the exponentiated saddle point action (31), leads to the result (6) scaled by the regularization-dependent prefactor (γt H ) −1 . The latter may be considered as an artifact of the path integral normalization and drops out in all physical observables normalized by the partition sum.

Partition function from Liouville quantum mechanics
We next discuss different approach to computing observables which is tailored to handle regimes of large fluctuations. The idea is to map the Liouville path integral to an equivalent Schrödinger equation and solve the latter. Where the partition function is concerned, this procedure leads to results identical to those of the stationary phase approach. The exactness of the latter in connection with Z has indeed been stated before [10], although we are not aware of a proof and the partition sum does not appear to satisfy standard criteria [28] for the semiclassical exactness of path integrals in an obvious ways. At any rate, pronounced deviations between the two approaches will be observed when other quantities are considered.
The functional integral in Eq. (24) lends itself to a straightforward quantum mechanical interpretation. To this end, consider the Liouville Hamiltonian operator and the corresponding imaginary time propagator: Comparison with Eq. (24) shows that the regularized Schwarzian partition function can be identified with the quantum propagator provided we have the identification We here require that at the integration boundaries, φ(±β/2) φ (±β/2) the semiclassical approximation remains valid, including in regimes where fluctuations are otherwise large. The justification for this assertion is that in the limit t H → ∞ the boundary field amplitudes assume large values and fluctuations are relatively less pronounced than in the interior of the integration interval. Under these conditions, the semiclassical expressions (28), (29) yield Eq. (35), which in turn leads to the equivalence Π ∼ Z.
The HamiltonianĤ, Eq. (32), has a continuous spectrum labeled by a quantum number k which may be interpreted as the momentum conjugate to φ. Its energy eigenvalues are given by E k = k 2 /2M, and are independent of γ. The dependence on this parameter is in the normalized wave functions which take the form where K α are modified Bessel functions, and Γ is the Gamma function. With these functions the spectral representation of the propagator (33) reads as Into this expression we substitute φ 1 = φ 2 = φ 0 , which means that we have to consider the (real) wave function amplitudes where the large argument asymptotic of the Bessel functions was used. Noting that (Γ(2ik)Γ(−2ik)) −1 = 2k sinh(2πk)/π the partition function (34) assumes the form The last expression implies that the average many-body DoS, ρ(E), of the SYK model is given by (7). Here, E measures the positive energy of reparameterization fluctuations and ρ(E) is the corresponding density of states relative to the ground state energy. This remarkably simple expression was obtained in Ref. [10] in the limit of a large number of interacting fermions and more recently in Ref. [13] by the combinatorial analysis of averaged moments of the Hamiltonian operator. We here show how it follows from the Liouville partition function. The square-root singularity of the DoS at low energies resembles the behavior of effective random matrix models close to their mean field spectral band gaps. In the case of SYK the corresponding square root behavior is observed at the energy scale ∼ 1/M ∼ J/N log N, where the reparameterization fluctuations become strong. As we explain in the Introduction it is this branch cut singularity of the many-body DoS which is responsible for the universal power-law behavior of the correlation functions at long times.
Finally, performing integration over k in Eq. (39), we obtain the partition sum (6) in exact agreement with the Gaussian approximation.

Out of time correlation function I: path integral representation
In this and the following sections we discuss the dynamics of the OTO correlation function. Unlike with the partition sum, it will turn out that now the proper treatment of fluctuations beyond the Gaussian level becomes essential. Let us consider the general four point function where the angular brackets denote both, disorder averaging, and the quantum thermal average tr(exp(−βĤ)(. . . )), and T τ is imaginary time ordering. From this function, the OTO correlation function (3) is obtained by the specific choice of complex time arguments (see Fig. 3) at s = 0, and the time ordering is along the contour [29] shown in Fig. 3. Here we used the fact that the correlation function depends only on differences between neighboring times to switch to a symmetric arrangements of the arguments ±it/2 relative to the real axis. (This configuration turns out to be more convenient in the subsequent integration over fluctuations.) Throughout, we will largely work with real time arguments, τ i ∈ [−β/2, β/2] and perform the analytic continuation to finite imaginary increments ±it/2 only at the final stages of the calculations. The shift parameter, s, in Eq. (41) enters the stage when the position of the reparameterization singularity at τ * , Fig. 2, relative to the observation times becomes of importance. In such cases, we shift τ * → −β/2 (which is always possible on account of the periodicity of the imaginary time theory) and offset the observation times by a parameter, s, which then is integrated over. We will work under the assumption that only soft reparameterization fluctuations are essential to the behavior of this function. The four point function then assumes the form of a product of two two point functions, where G β ([ f ]) is the reparameterized Green function (17) and the functional average defined by Eq. (24) leads to correlations between them. To bring the Green functions into a form suitable for functional integration, we substitute (21) into (17) and use (12) to obtain This representation will be an important building block in all subsequent calculations. For later reference we note that the integration over the auxiliary variable comes with a convergence generating factor, ∼ exp(−αδ), where δ ∼ J −1 is of the order of the inverse UV cutoff. The reason is that we are operating within the framework of an effective low energy theory and times of the order ∼ J −1 cannot be effectively resolved. This translates to a smearing of the same order in the arguments of the Green functions, h(τ 1 ) − h(τ 2 ) + δ, which after exponentiation acts as a convergence generator. In this form, the Green functions can now be substituted into Eq. (24) and we obtain the representation where and S [φ] is given in Eq. (33). Eq. (44) and (26) define the path integral we need to compute. Its exponent suggests an interpretation in terms of an effective quench dynamics: 2 between the times τ 4 , τ 2 ,τ 3 and τ 1 , respectively, the strength of the Liouville potential effectively jumps as γ, γ + α 2 , γ + α 1 + α 2 , γ + α 1 , γ (cf. the right panel in Fig. 3). In the following sections we analyze this kink dynamics by stationary phase methods (short times/high temperatures) and via the corresponding Schrödinger equation (long times/low temperatures), respectively.

OTO Correlation function II: short times/high temperatures
The computation of the OTO path integral by stationary phase methods parallels that employed in section 4 for the partition sum. To simplify the subsequent calculations we again set β ≡ π. We also re-introduce the scaled integration variable (25) and in addition scale the auxiliary integration variables as This leaves the integral invariant except that the action now assumes the form (26).
Our later analysis will self-consistently show that the characteristic values of the integration variables, α i = O(1) MT . This means that 'quench potentials' α i e ϕ are weak as compared to the unperturbed Liouville potential e ϕ , and that it makes sense to expand around a stationary phase solution δ ϕ S [ϕ] = M(−ϕ + 2e ϕ ) = 0 ignoring the former. The solution to this equation is given byφ = − ln(cos 2 (ϕ)), i.e. by Eq (28) neglecting the shift parameter λ which is inessential except for an infinitesimal neighborhood of the integration boundaries. If we evaluate the functional integral on this configuration, the integration over the auxiliary variables yields, where we defined and in the final expression encounter the uncorrelated product of two mean-field Green functions (15). This suggests that fluctuations around the mean field will have two distinct effects: they will modify the form of the individual Green functions and, more importantly, generate correlations between them. To disentangle these influences, we introduce the ratio where G(τ, τ ) are the mean field Green functions corrected by fluctuations. The leading (in M −1 ) contribution to the fluctuation action comes from the quadratic expansion of the action, S [φ + δϕ] in δϕ, In addition, the fluctuation offset δϕ appears in the pre-exponential factors, and in the quench potentials α i e ϕ . These multiple appearances make the full integration over fluctuations technically cumbersome. It turns out, however, that for large M fluctuations of the pre-exponential factors exp(δϕ(τ i )/4) contribute only negligibly to the both the renormalization of the Green functions and to correlations, and may be neglected. The expansion of the quench actions where terms with superscript (n) are of nth order in δϕ, produces two contributions, S (2) α i − 1 2 (S (1) α i ) 2 , where the angular brackets denote functional averaging over the fluctuation action (50). These lead to a weak renormalization of the ith Green function G β which drops out when the ratio f is built. The term of interest which does generate non-vanishing correlations is the functional average of the cross ratio The straightforward but lengthy computation of the functional expectation value is detailed in Appendix D and yields I(τ 2 , τ 4 ) − sin(4τ 2 ) cos(2τ 4 ) cos 2 τ 4 2M sin 2 (2τ 2 ) sin 2 (2τ 4 ) where we used that (including during the analytic continuation) we may impose the constraint τ 3 = τ 4 + π/2 and τ 1 = τ 2 + π/2, and the sign indicates the omission of terms exponentially vanishing under under the continuation to large imaginary times it iβ. Using this result, we find that the leading contribution in M −1 to the correlation ratio Eq. (49) is given by Substituting Eq. (48), the α-integrals appearing in this expression become and analogously for the second integral. We substitute them in Eq. (52) to obtain where in the final step we have performed the analytic continuation to the time arguments (cf. Eq. (41)) τ 2 → −π/4 − it/2, τ 4 → −π/2 + it/2. Expressed in the original language of correlation functions, and re-introducing the temperature parameter, this result assumes the form This agrees, including pre-factors, with Eq. (6.59) of Ref. [6], where it was obtained from the Schwarzian action, without recourse to its Liouvillian reformulation. This expression may be trusted up until the second term in the brackets becomes comparable with unity, which happens at the Ehrenfest time t E , given by Eq. (4).

OTO correlation function III 4: intermediate times/high temperatures
In this section we briefly discuss what happens for times larger than the Ehrenfest time t E , Eq. (4), at which the correction discussed previously overpowers the leading order result (56). Yet we restrict ourselves to times smaller than ∼ M, where the strong Goldstone mode fluctuation regime is entered. This intermediate regime is still amenable to stationary phase methods. However, the technicalities get a little complicated because the quench potentials α i e ϕ can no longer be neglected in the solution of the stationary phase equations. Referring to Appendix F for details, the stationary phase procedure leads to the intermediate result for OTO correlation with fixed α 1,2 : .
This formula already contains the essence of the result: exponential decay of the correlation function with an exponent whose real part scales as ∼ exp(−t) → exp(−tπ/β). However, if we aim to fix the result including pre-exponential factors, the integration over auxiliary variables needs to be performed. Here, we recalled that the α-integrations come with a convergence generating factor δ ∼ J −1 (cf. discussion after Eq. (43)) which under the rescaling by γ ∼ J gets promoted to a factor of O(1). This factor acts to regularize an otherwise logarithmically divergent integral. The integration is not entirely trivial. One observes that the integrand does not couple to physical parameters, except through the variable combination ω(α 1 , α 2 ). This suggests to introduce a new variable, λ = λ(α 1 , α 2 ) ≡ ln((1 + α 1 )(1 + α 2 )/(1 + α 1 + α 2 )). Expressed in terms of λ, the exponent becomes Gaussian, and the rest of the integral, dα 1 dα 2 (. . . )δ(λ(α 1 , α 2 )−λ) ≡ ρ(λ) turns into an effective 'spectral density'. Using that for large M only small deviations of λ off zero are relevant, one may verify that ρ(λ) ln(λ −1/2 )/ √ λ, and the integral becomes Finally, re-introducing β and recalling the value of the finite temperature Green functions G β (β/2) ∼ (βJ) −1/2 , we obtain where the definition of the Ehrenfest time, exp(πt E /β) = M/β was used. To exponential accuracy this agrees with the black hole shock wave analysis of Ref. [6] which predicted a decay ∼ t exp(−πt/β). For completeness we mention that the above discussion assumed a fixation of the real parts Re(τ i ) of the time arguments relative to that, τ * = β/2, of the singularity which is necessarily present in a bijective transformation [−β/2, β/2] → R (cf. discussion in the beginning of section 3.) This fixation may be abandoned by introducing a shift parameter τ i → τ i + s and integrating over the latter (cf. Fig. 3.) On top of that, one may place the singularity in-between any of the arguments τ i , and sum over all four choices. While the relative ordering of the singularity and the observation times turns out to be irrelevant, the continuous shift parameter begins to play a role when we turn to the complementary regime of large observation times:

OTO correlation function IV: long times/low temperatures
The previous sections served benchmarking purposes where we reproduced known results within the framework of the finite temperature extension of Liouville theory. We now turn to uncharted territory and address the theory at low temperatures, where quantum fluctuations are large and different methodology is required. As in section 5, the strategy is to turn a foe into a friend and use that largeness of the φ-fluctuations implies confinement of the fluctuations of its conjugate momentum, k.
Once more the calculation simplifies if a convenient set of complex time arguments is chosen. Presently, we find it advantageous to start from the configuration (41) at it = τ ∈ R. Here, τ serves as a real time parameter which will eventually be continued into the complex plane, and s ∈ [0, β/4] is the parameter controlling the relative position of the observation times and the time, τ * , of the reparameterization singularity (cf. discussion at the end of the previous section). The correlation function is then obtained as where G 4 is given by Eq. (44), and the integral implements an average over the shift parameter, s. The ellipses stand for the three contribution of different ordering of the singularity relative to the observation times. One can explicitly verify that all four orderings give identical result for the real part while a pairwise cancellation of imaginary parts occurs. Therefore we will focus on the one explicitly displayed throughout an account for the remaining ones by a multiplicative factor of 4. It will be convenient to scale the integration variable in the path integral as φ → φ + ln γ, and the auxiliary variables as α i → α i /γ. This operation leaves the form of the path integral invariant but changes the strength of the potentials in the Liouville action as γ → 1, and α i → α i /γ. As in section 5, we re-interpret the path integral in a Schrödinger picture where it describes time evolution under the Hamiltonian Eq. (32). However, the presence of the quench potentials α 1,2 e ϕ implies that the strength of the Liouville potential jumps from 1 to 1 + α 1 /γ and 1 + α 2 /γ in the time intervals [τ 2 , τ 1 ] and [τ 4 , τ 3 ], respectively (cf. Eq. (45).) This requires the insertion of different sets of eigenfunctions labeled as |k , |k (α) , depending on whether the quench potentials are absent or present with value α, respectively.
Specifically, the presence of four observation times and one singularity requires the insertion of five resolutions of unity. The spectral representation of the correlation function then assumes the form where τ i j = τ i − τ j are the time differences Central to the expression above are the matrix elements, k (α) |e φ 4 |k (α ) between eigenstates of different Liouville potential. In Appendix E we show that for the small momenta k relevant to the correlation function at large times they show a high degree of universality. The matrix elements factorize into products of the α-independent normalization factors N k = 2 Γ(2ik) (cf. Eq. (36)) and M −1/4 times a term which depends only on α. The integral over these variables can thus be performed to yield a numerical factor C = O(1). This leads to the intermediate result Here, we used that for the diverging boundary conditions, φ 0 → ∞, the boundary wave functions, , converge to a product of the normalization factor and a k-independent term. The latter cancels against the corresponding factor appearing in the partition sum (see below) and is ignored.
Using the expansion |N(k)| 2 = |Γ(2ik)| −2 4k 2 , and the auxiliary relation dk 2π k 2 exp(−E k τ) = 1 √ 2π (M/τ) 3/2 , we obtain the proportionality, Here, the ellipses · · · ∼ M 0 ds(. . . ) refer to the contribution to the integral where the observation points are close to the boundary. The short time distance Green functions appearing in this regime can no longer be described within the low momentum approximation used above. However, referring to Ref. [21] for details we note that they behave as ∼ 1/Ms 1/2 , matching the above s −3/2 law at s ∼ M. In this way the boundary contributions effectively regularize the s −3/2 singularity of the integral. We also note that the accumulated appearance of τ −3/2 scaling in this relation is a hallmark of Liouville quantum mechanics. As mentioned in the Introduction, cf. Eq. (8), it may be qualitatively understood assuming statistical independence of many-body energies and matrix elements and utilizing square root singularity in the many-body DoS.
Finally, noting the scaling of the partition sum Z(β) ∼ dk| φ 0 |k | 2 exp(−βk 2 /2M) ∼ (M/β) 3/2 , substituting the time arguments, Eq. (63), and continuing to real times, τ → it, we obtain Incidentally, the t −6 power law scaling coincides with that found for non-interacting spinless fermions [31]. The application of the same spectral decomposition procedure to the Greens function, G(τ), yields the result Using this formula to normalize the OTO function, F(t) → F(t)/[G(β/2)] 2 , we obtain the result describing regime 4 in Eq. (5). The results derived above can readily be extended to the regime of long times but high temperatures, β M t. To this end we substitute into Eq. (65) the high temperature expression for the partition sum (6), and use the scaling G(β/2) ∼ (Jβ) −1/2 of the 2-point Greens function to normalize. As a result, Eq. (66) changes to which now describes the OTO function in the regime 3 of Eq. (5).

Discussion
In this paper we explored the role played by large conformal Goldstone mode fluctuations of the SYK model within a finite temperature framework suitable for the study of out of time order correlations. Our theory is constructed by exploiting the extensive time-reparameterization freedom of the SYK model. This allows to describe Goldstone mode fluctuations in terms of a real valued field defined on the periodic imaginary time interval which, however, necessarily contains a point of singularity. The handling of this singularity required some care and we tested the theoretical framework by comparison to various known, or conjectured results, including the temperature dependence of the many-body partition sum, the exponential buildup of OTO decorrelations at short times, and their exponential decay at times larger than the Ehrenfest time.
The truly novel content of the theory unfolds at times larger than a scale M ∼ N ln(N)J −1 , where N is the number of particles in the system or, equivalently, a logarithmic measure for the dimension ∼ 2 N/2 of the many-body Hilbert space. This regime is inaccessible to theoretical approaches taking thermodynamic limits N → ∞ at a fixed time, but arguably may play a key role in the description of the regularizing effects of quantum mechanics (via finite N) on long time fluctuations. Our main finding is that in such regimes, and equally in regimes of low temperatures, T < M −1 at generic times t > T −1 , four-point OTO correlations decay with the universal power law ∼ t −6 .
While the calculations required to nail this behavior are technical in parts (mainly due to the required careful handling of the singularity), the 6 = 4 × 3/2 power law has its origin in the celebrated t −3/2 universality of the Liouville quantum mechanics [26]. It is a robust feature of this system that temporal correlations of local operators decay as Ô (0)Ô(t) ∼ t −3/2 . This power-law may also be interpreted in terms of the √ E scaling of the average low-energy many-body DoS, Eq. (7), as discussed in the Introduction. The quadrupling of the 3/2 exponent has to do with the definition of the OTO correlation function whose particular time ordering of operators requires the introduction of four temporal contours, Fig. 3. This time arrangement makes all four time arguments of the correlation function (40), (41) to be separated by long intervals ∼ |t ± iβ/4|. One may generalize this observation and claim that any multipoint OTO correlation is given by the product of 3/2 powers of all long ( M) time intervals involved.
Note added: Recent preprint [32] discusses one-loop (semiclassical) exactness of the partition sum of the Schwarzian theory (6) using the Duistermaat-Heckman formula [33]. It remains to be seen if this technique is also helpful in evaluation of correlations functions. lim n→0 Z n = 1. Averaging over disorder, we obtain To pass from the abbreviation,G, for a Grassmann bilinear to an integration variable, G, we introduce in the functional integral, which entitles us to replaceG → G in the quartic term. As a result, the action has become quadratic in Grassmann variables, and the Gaussian integration over these produces the Pfaffian ∂ τ δ ab + Σ ab τ,τ The situation simplifies further if we observe that eventually an analytic continuation to large imaginary values it will be carried out. This implies that in the integration over trigonometric functions only those contributions need to be kept that will not vanish exponentially. Under this condition the straightforward but tedious calculation of the auxiliary integral simplifies to the result (52).
Appendix E. Computation of the matrix elements k (α) |e φ 4 |k (α ) After the scaling, γ → 1, α → α/γ the Liouville eigenfunctions in the φ-representation are given by Eq. (36) with the replacement γ → 1 + α. As a consequence, these matrix elements of the operator exp(φ/4) assume the form For generic values of k, k the integral evaluates to an complicated configuration of hypergeometric functions. However, at large times only small momenta matter, and to zeroth order in k the result simplifies to W(k, k |α, α ) N k N k M 1/4 W(α, α ), W(α, α ) = Γ 2 ( 1 4 ) 2 5/4 where K is the complete elliptic integral of the first kind. For small values of α the functions K are regular, and for large arguments weakly decay as α −1/4 . This translates to a logarithmic singularity of the α-double integral. What comes to rescue is the convergence generating factor (cf. discussion below Eq. (43)). In the scaled framework, γ → 1, this factor is upgraded as δ ∼ J −1 to δγ ∼ 1 and to the integral we need to perform reads ∞ 0 dα 1 dα 2 √ α 1 α 2 W(0, α 2 )W(α 2 , α 1 + α 2 )W(α 1 + α 2 , α 1 )W(α 1 , 0)e −α 1 −α 2 ≡ C, where C is a constant of order unity. auxiliary variables as α i → 2Mα i . The entire action then is multiplied by a factor M, and that the functional integral picks up the same multiplicative factor. The scaling of variables also implies that the functional integral for the four point function picks up global multiplicative factor 2M. Due to the piecewise constancy of the potential strength, the solutions to the stationary phase equationsφ − 2γ i exp(φ) = 0, γ 0 = 1, γ 1 = 1 + α 1 , γ 2 = 1 + α 1 + α 2 , γ 3 = 1 + α 2 , γ 4 = 1 for the five temporal regimes defined by the time arguments −β/2, τ 4 , τ 2 , τ 3 , τ 1 , β/2 assume the form ϕ = ln(ω 2 γ −1 i cos −2 (ωτ + δ)), where ω and δ are free parameters. Requiring the usual boundary singularityφ(±π/2) → +∞ and continuity of the solution in the bulk a straightforward fixation procedure yields for these parameters δ 0 = −δ 4 = i 2 ln 1 + α 1 + α 2 (1 + α 1 )(1 + α 2 ) , where ω is specified in Eq. (57). (As a technical remark we note that the evaluation of the matching conditions is facilitated by using that after analytic continuation, the matching points have large imaginary parts, |Im(τ i )| = t/2, and an approximation ϕ = ln(ω 2 γ −1 i ) ± 2i(ωτ − δ) is possible, where the sign is postive/negative for τ 4,3 /τ 1,2 .) The substitution of these configurations into the pre-exponential factor yields where Eq. (41) was used. Similarly, the substitution of the stationary solution into the action yields S [φ] = −2Mπω 2 0 . (To derive this result, the regularization of the solution near the boundary by the parameter λ of Eq. (28) needs to be taken into account. As before it leads to a cancellation of the formally divergent factor γt H in the action.) Finally, we note that the integration of fluctuations around the optimal configuration leads to a fluctuation determinant independent of the observation times. This factor cancels against the identical pre-exponential factor multiplying the partition sum (6) and may be ignored. Combining terms, we arrive at Eq. (57).