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

Introduction

The Sachdev-Ye-Kitaev (SYK) model is a system of $N$ (Majorana) fermions, $\chi_{i}$ with $i=1,\ldots,N$, subject to a four-fermion interaction

with Gaussian distributed random matrix elements $J_{ijkl}$ with zero mean and a variance given by $\langle|J_{ijkl}|^{2}\rangle=6J^{2}/N^{3}$. The seeming simplicity of this model is deceptive. At low excitation energies it exhibits an asymptotically exact conformal symmetry which manifests itself in the infinite-dimensional freedom to re-parameterize time, $t\to f(t)$, in the description of long time correlations. This ‘nearly conformal symmetry’ (NCFT) 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 .

In Ref. 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^{\prime}(t)=\exp(\phi(t))$ whereupon the Goldstone mode integral assumes the form of a path integral over $\phi$, with an (approximately) time local action coinciding with that of so-called Liouville quantum mechanics . Standard methods of quantum mechanics are then applicable to perform the integration, including the regime where the $\phi$-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\gg 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 , $G(\epsilon)\sim|\epsilon|^{-1/2}$ has a divergent amplitude at low energies, the inclusion of Goldstone mode fluctuations shows that $G(\epsilon)\sim|\epsilon|^{1/2}$ for energies $|\epsilon|\ll 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. we analyzed the effect of Goldstone modes within a zero temperature framework. However, the majority of observables describing the dynamical behavior of the SYK-model — four point functions in general, and out-of-time-order (OTO) correlation functions in particular variants of the partition sum , 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 , $\sim\exp(2\pi Tt)$, 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

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\gg M^{-1}$, we need to discriminate between short and large times, $t\ll M$ and $t\gg M$, respectively. In the short time regime, Goldstone mode fluctuations are weak, and one identifies two different regimes characterized by the exponential loss of correlations, depending on whether $t<t_{E}$, or $t>t_{E}$, where

plays the role of an Ehrenfest time in the problem. However, at large times, $t\gg M$, fluctuations are strong and generate universal power law scaling, $F(t)\sim t^{-6}$. For small temperatures, $T\ll M^{-1}$, Goldstone mode fluctuations affect the picture throughout the entire domain, no exponential regimes are found, and the power law, $F(t)\sim t^{-6}$, holds for all times $t\gg 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 to the regime of long times/low temperatures. These correlations appear as a robust consequence of 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\geq 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 .

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\ll J$ scales as

in agreement with the results of Ref. . This formula is directly related to the many-body density of states (DoS), which we find for energies $E\ll 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\sim\sqrt{N}$, and by statistical analysis of moments of the random interaction operator . In our present analysis, the DoS reflects the energy stored in Goldstone mode fluctuations. We also note that at high temperatures, $T\gg 1/M$, or energies, $E\gg 1/M$, the constant $M$ in Eqs. (12) and (13) acquires a week logarithmic dependence on either $T$ or $E$, respectively (see more comments on this in section 3).

Interestingly, the result (12) for the partition sum entails an heuristic explanation for the power-law formation in the OTO correlation function (10) 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\rangle$ in (3) to obtain the ‘Lehmann representation’

Now, let us assume that the Majorana matrix elements $\langle m|\chi_{i,j}|n\rangle$ and the eigenenergies, $\epsilon_{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 may contain some truth. If so, and if at the characteristic energy differences $|\epsilon_{m}-\epsilon_{m}|\sim t^{-1}\sim M^{-1}$ much larger than the many-body level spacing $\sim e^{-N/2}$ statistical correlations between the levels of different spectral sums can be neglected, the correlation function simplifies to

where Z(s)=\big{\langle}\sum_{n}e^{-s\epsilon_{n}}\big{\rangle}=\int dE\,\rho(E)\,e^{-sE} with ${\rm Re}(s)>0$ is the Laplace transform of the average many-body DoS, Eq. (13), which at $s=\beta$ 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 $\rho(E)\sim E^{1/2}$, it leads to $Z(s)\sim s^{-3/2}$ at $|s|\gg M$. This is exactly the 3/2 law, observed in the long-time behavior of two-point functions . It in turn implies $F(t)\sim t^{-6}$ as announced in Eq. (10).

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. (10) 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 after the averaging over disorder and integration over Grassmann-fermion degrees of freedom have been performed (see A)

Here, $G^{ab}_{\tau,\tau^{\prime}}$ and $\Sigma^{ab}_{\tau,\tau^{\prime}}$ are matrix fields carrying an index structure, $a,b=1,\dots,n$ in a space of $n$ replicas. The first term of the action resembles the free energy of a fermion system and identifies $\Sigma$ as a dynamical self energy. The second term reflects that the average over disorder, $\langle\hat{H}^{2}\rangle$ 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 $\Sigma,G$.

The action (16) exhibits an exceptionally high level of symmetry, provided $\partial_{\tau}$ term is neglected: it is then invariant under reparameterizations of time, $\tau\to f(\tau)$, where $f(\tau)$ may be any invertible and differentiable function. (Invertibility requires globally monotonicity and without loss of generality we assume $f^{\prime}(\tau)>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 $\Sigma$ become more tangible at the mean-field level, where the stationary phase equations 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}\equiv\big{[}G^{ab}_{\tau,\tau^{\prime}}\big{]}^{3}. These equations can be solved in the long time limit $\tau-\tau^{\prime}\gg 1/J$, where $\partial_{\tau}$ may be neglected. A configuration, $G_{\infty},\Sigma_{\infty}$ solving the equations at zero temperature, $\beta\to\infty$, reads as

where $b=(4\pi)^{-1/4}$. However, this solution is not unique . Under the action of the symmetry operation it transforms as

The transformation (19) 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 $[-\beta/2,\beta/2]$ : the reparameterization

provides an invertible map from the open imaginary time interval $[-\beta/2,\beta/2]$ onto the entire axis. According to Eq. (19) the corresponding Green function reads

and similarly for $\Sigma_{\beta}(\tau_{1}-\tau_{2})$. This defines a periodic and time translationally invariant solution of the finite temperature saddle-point equation (17). The group property of the symmetry implies that reparameterization of compactified time $\tau\to f(\tau)$, mapping the interval $[-\beta/2,\beta/2]$ onto itself and obeying the periodicity constraint

generate a family of finite temperature solutions $G_{\beta}([f],\tau,\tau^{\prime})$, periodic on the time interval but lacking translational invariance (i.e. the functions $G_{\beta}[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 soft-mode manifold parameterized by the functions $f$ identified in the previous section. In Ref. 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.

Recalling that $G_{\beta}=G_{\infty}([g])$, we have $G_{\beta}[f])=G_{\infty}([g\circ f])=G_{\infty}[h]$, or

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 $\tau,\tau^{\prime}\in[-\beta/2,\beta/2]$.

Likewise, any finite-temperature observable, ${\cal O}_{\beta}[f]$, formulated in terms of reparameterizations, $f$, of the finite temperature soft mode manifold may be transformed into the zero temperature observable ${\cal O}_{\infty}[h]$, parameterized by $h=g\circ f$. Expectation values of such observables are obtained by integration over reparameterizations,

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 $\Delta$ defined self-consistently via $\Delta=1/M(\Delta)$, which with log accuracy is resolved as $\Delta=64\sqrt{\pi}J/(N\ln N)$. At smaller energies, $E<\Delta$, the Goldstone mode fluctuations proliferate and $M(E)$ saturates to the constant value $M\equiv 1/\Delta$ 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={\rm min}(T,1/t)$, where $t$ is a time variable in the correlation function.

Before turning to the actual integration procedure in Eq. (25), we impose one more change of variables and introduce the degree of freedom, $\phi(\tau)$, of Liouville quantum mechanics as

where the positivity $h^{\prime}(\tau)>0$ guarantees that $\phi(\tau)$ is real. The key advantage of this parameterization is the flatness of its integration measure, $\mu[h]\,{\cal D}h={\cal D}\phi$ . However, attention must be payed to the singularity of all admissible $h$-functions at $\tau=\pm\beta/2$. It translates into a diverging integral $\int_{-\beta/2}^{\beta/2}d\tau\,h^{\prime}(\tau)=\int_{-\beta/2}^{\beta/2}d\tau\,e^{\phi(\tau)}=\infty$. 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}\to\infty$ at the end of the calculation. Trajectories $\phi(\tau)$ satisfying this condition thus become regularized at the boundaries of the time interval, $\phi(\pm\beta/2)=\phi_{0}$ with $\phi_{0}\to+\infty$ as $t_{H}\to\infty$.

In the $\phi$-language the soft mode action reads as

and expectation values of observables ${\cal O}_{\beta}$ assume the form

where $\gamma$ is a Lagrange multiplier enforcing the condition (29), and $Z_{\beta}=\int\mathcal{D}\phi\exp(-S[\phi])$ is the partition sum. Notice that the full derivative $\phi^{\prime\prime}$ cannot be ignored since, due to singularity discussed above, $\phi^{\prime}(\pm\beta/2)\to\pm\infty$. It is this term which neutralizes the formally infinite constant $\gamma t_{H}$. The same singularity also shows in the boundary conditions $\phi(\pm\beta/2)=\phi_{0}$. Our treatment below establishes a balance between the parameters $t_{H}$, $\phi_{0}$, and $\gamma$ such that the final answers do not depend on regularization specifics.

Saddle point calculation

At sufficiently high temperature, and/or sufficiently short time intervals $\tau-\tau^{\prime}$ the functional integral (31) 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 $\beta\equiv\pi$. In the final results the $\beta$-dependence can then be re-introduced by scaling all quantities with dimension of time as $\tau\to\tau\pi/\beta$. Expressed in the new language, the action assumes the form

and the equation of motion reads as $\varphi^{\prime\prime}(\tau)=2e^{\varphi(\tau)}$. In the limit $t_{H}\to+\infty$ this is solved by $\bar{\varphi}_{\infty}(\tau)=-\ln\left(\cos^{2}\tau\right)$. (The dependence on $t_{H}$ is implicit through Eq. (34), which for $t_{H}\to\infty$ implies diverging boundary conditions $\bar{\varphi}_{\infty}(\tau\to\pm\pi/2)\to\infty$.) The meaning of this solution becomes apparent when we formulate it in the language of the $h$-representation, $\bar{h}_{0}(\tau)=\frac{2M}{\gamma}\int d\tau e^{\bar{\varphi}_{\infty}(\tau)}=\frac{2M}{\gamma}g(\tau)$, where $g$ is the $\tan$-function of Eq. (21). It is straightforward to check that for this function $G_{\infty}([\bar{h}_{0}],\tau,\tau^{\prime})=G_{\beta}(\tau-\tau^{\prime})$ reduces to the familiar translationally invariant finite temperature solution of the Dyson equation (22).

Finite values of $t_{H}$ regularize the infinity in the boundary conditions of $\bar{\varphi}_{\infty}$. This regularization can be achieved by shifting the argument of the $\cos$ in the solution as,

For small $\lambda$ we find that the integral in (34) evaluates to $4/(\pi\lambda)$, and this fixes the shift parameter as

For this configuration, we find a boundary value, $\exp(\bar{\varphi}(\pm\pi/2))\equiv\exp(\varphi_{0})=(2/\pi\lambda)^{2}$.

Finally, a straightforward calculation shows that the stationary phase action is given by

Gaussian fluctuations around $\bar{\varphi}$ lead to a fluctuation determinant whose calculation is detailed in C. Multiplying this factor with the exponentiated saddle point action (38), leads to the result (12) scaled by the regularization-dependent prefactor $(\gamma 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 , although we are not aware of a proof and the partition sum does not appear to satisfy standard criteria 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. (31) 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. (31) shows that the regularized Schwarzian partition function can be identified with the quantum propagator

We here require that at the integration boundaries, $\phi(\pm\beta/2)\simeq\bar{\phi}(\pm\beta/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}\to\infty$ 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 (35), (36) yield Eq. (42), which in turn leads to the equivalence $\Pi\sim Z$.

The Hamiltonian $\hat{H}$, Eq. (39), has a continuous spectrum labeled by a quantum number $k$ which may be interpreted as the momentum conjugate to $\phi$. Its energy eigenvalues are given by $E_{k}=k^{2}/2M$, and are independent of $\gamma$. The dependence on this parameter is in the normalized wave functions which take the form

where $K_{\alpha}$ are modified Bessel functions, and $\Gamma$ is the Gamma function. With these functions the spectral representation of the propagator (40) reads as

Into this expression we substitute $\phi_{1}=\phi_{2}=\phi_{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 $(\Gamma(2ik)\Gamma(-2ik))^{-1}=2k\sinh(2\pi k)/\pi$ the partition function (41) assumes the form

The last expression implies that the average many-body DoS, $\rho(E)$, of the SYK model is given by (13). Here, $E$ measures the positive energy of reparameterization fluctuations and $\rho(E)$ is the corresponding density of states relative to the ground state energy. This remarkably simple expression was obtained in Ref. in the limit of a large number of interacting fermions and more recently in Ref. 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 $\sim 1/M\sim 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. (46), we obtain the partition sum (12) 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

at $s=0$, and the time ordering is along the contour 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 $\pm 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, $\tau_{i}\in[-\beta/2,\beta/2]$ and perform the analytic continuation to finite imaginary increments $\pm it/2$ only at the final stages of the calculations. The shift parameter, $s$, in Eq. (48) enters the stage when the position of the reparameterization singularity at $\tau^{\ast}$, Fig. 2, relative to the observation times becomes of importance. In such cases, we shift $\tau^{\ast}\to-\beta/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_{\beta}([f])$ is the reparameterized Green function (24) and the functional average defined by Eq. (31) leads to correlations between them. To bring the Green functions into a form suitable for functional integration, we substitute (28) into (24) and use (18) 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, $\sim\exp(-\alpha\delta)$, where $\delta\sim 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 $\sim J^{-1}$ cannot be effectively resolved. This translates to a smearing of the same order in the arguments of the Green functions, $h(\tau_{1})-h(\tau_{2})+\delta$, which after exponentiation acts as a convergence generator.

In this form, the Green functions can now be substituted into Eq. (31) and we obtain the representation

Eq. (52) and (33) define the path integral we need to compute. Its exponent suggests an interpretation in terms of an effective quench dynamics: Quantum quench is a protocol implying that parameters of a system are suddenly changed, thus at times $t<0$ a unitary evolution of the systems proceeds according to the Hamiltonian $\hat{H}_{0}$ while at times $t>0$ it is governed by the perturbed Hamiltonian $\hat{H}\neq\hat{H}_{0}$. In general, such sudden change of $\hat{H}_{0}$ is not supposed to be weak. See e.g. Ref. for the recent review. between the times $\tau_{4}$, $\tau_{2}$,$\tau_{3}$ and $\tau_{1}$, respectively, the strength of the Liouville potential effectively jumps as $\gamma,\gamma+\alpha_{2},\gamma+\alpha_{1}+\alpha_{2},\gamma+\alpha_{1},\gamma$ (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 $\beta\equiv\pi$. We also re-introduce the scaled integration variable (32) and in addition scale the auxiliary integration variables as

This leaves the integral invariant except that the action now assumes the form (33).

Our later analysis will self-consistently show that the characteristic values of the integration variables, $\alpha_{i}=\mathcal{O}(1)\ll MT$. This means that ’quench potentials’ $\alpha_{i}e^{\varphi}$ are weak as compared to the unperturbed Liouville potential $e^{\varphi}$, and that it makes sense to expand around a stationary phase solution $\delta_{\varphi}S[\varphi]=M(-\varphi^{\prime\prime}+2e^{\varphi})=0$ ignoring the former. The solution to this equation is given by $\bar{\varphi}=-\ln(\cos^{2}(\varphi))$, i.e. by Eq (35) neglecting the shift parameter $\lambda$ 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,

and in the final expression encounter the uncorrelated product of two mean-field Green functions (22). 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(\tau,\tau^{\prime})$ 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[\bar{\varphi}+\delta\varphi]$ in $\delta\varphi$,

The straightforward but lengthy computation of the functional expectation value is detailed in D and yields

where we used that (including during the analytic continuation) we may impose the constraint $\tau_{3}=\tau_{4}+\pi/2$ and $\tau_{1}=\tau_{2}+\pi/2$, and the $\simeq$ sign indicates the omission of terms exponentially vanishing under under the continuation to large imaginary times $it\gg i\beta$. Using this result, we find that the leading contribution in $M^{-1}$ to the correlation ratio Eq. (57) is given by

Substituting Eq. (56), the $\alpha$-integrals appearing in this expression become

and analogously for the second integral. We substitute them in Eq. (60) to obtain

where in the final step we have performed the analytic continuation to the time arguments (cf. Eq. (48)) $\tau_{2}\to-\pi/4-it/2$, $\tau_{4}\to-\pi/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. , 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 (65). Yet we restrict ourselves to times smaller than $\sim 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 $\alpha_{i}e^{\varphi}$ can no longer be neglected in the solution of the stationary phase equations. Referring to F for details, the stationary phase procedure leads to the intermediate result for OTO correlation with fixed $\alpha_{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 $\sim\exp(-t)\to\exp(-t\pi/\beta)$. 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 $\alpha$-integrations come with a convergence generating factor $\delta\sim J^{-1}$ (cf. discussion after Eq. (50)) which under the rescaling by $\gamma\sim J$ gets promoted to a factor of $\mathcal{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 $\omega(\alpha_{1},\alpha_{2})$. This suggests to introduce a new variable, $\lambda=\lambda(\alpha_{1},\alpha_{2})\equiv\ln((1+\alpha_{1})(1+\alpha_{2})/(1+\alpha_{1}+\alpha_{2}))$. Expressed in terms of $\lambda$, the exponent becomes Gaussian, and the rest of the integral, $\int d\alpha_{1}d\alpha_{2}(\dots)\delta(\lambda(\alpha_{1},\alpha_{2})-\lambda)\equiv\rho(\lambda)$ turns into an effective ‘spectral density’. Using that for large $M$ only small deviations of $\lambda$ off zero are relevant, one may verify that $\rho(\lambda)\simeq\ln(\lambda^{-1/2})/\sqrt{\lambda}$, and the integral becomes

Finally, re-introducing $\beta$ and recalling the value of the finite temperature Green functions $G_{\beta}(\beta/2)\sim(\beta J)^{-1/2}$, we obtain

where the definition of the Ehrenfest time, $\exp(\pi t_{E}/\beta)=\sqrt{M/\beta}$ was used. To exponential accuracy this agrees with the black hole shock wave analysis of Ref. which predicted a decay $\sim t\exp(-\pi t/\beta)$.

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 $\phi$-fluctuations implies confinement of the fluctuations of its conjugate momentum, $k$.

where $G_{4}$ is given by Eq. (52), 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 $\phi\to\phi+\ln\gamma$, and the auxiliary variables as $\alpha_{i}\to\alpha_{i}/\gamma$. This operation leaves the form of the path integral invariant but changes the strength of the potentials in the Liouville action as $\gamma\to 1$, and $\alpha_{i}\to\alpha_{i}/\gamma$. As in section 5, we re-interpret the path integral in a Schrödinger picture where it describes time evolution under the Hamiltonian Eq. (39). However, the presence of the quench potentials $\alpha_{1,2}e^{\varphi}$ implies that the strength of the Liouville potential jumps from $1$ to $1+\alpha_{1}/\gamma$ and $1+\alpha_{2}/\gamma$ in the time intervals $[\tau_{2},\tau_{1}]$ and $[\tau_{4},\tau_{3}]$, respectively (cf. Eq. (53).) This requires the insertion of different sets of eigenfunctions labeled as $|k\rangle,|k^{(\alpha)}\rangle$, depending on whether the quench potentials are absent or present with value $\alpha$, 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 $\tau_{ij}=\tau_{i}-\tau_{j}$ are the time differences

Central to the expression above are the matrix elements, $\langle k^{(\alpha)}|e^{\frac{\phi}{4}}|k^{{}^{\prime}(\alpha^{\prime})}\rangle$ between eigenstates of different Liouville potential. In 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 $\alpha$-independent normalization factors ${\cal N}_{k}=\tfrac{2}{\Gamma(2ik)}$ (cf. Eq. (43)) and $M^{-1/4}$ times a term which depends only on $\alpha$. The integral over these variables can thus be performed to yield a numerical factor $C=\mathcal{O}(1)$. This leads to the intermediate result

Using the expansion $|\mathcal{N}(k)|^{2}=|\Gamma(2ik)|^{-2}\simeq 4k^{2}$, and the auxiliary relation $\int\frac{dk}{2\pi}\,k^{2}\exp(-E_{k}\tau)=\tfrac{1}{\sqrt{2\pi}}(M/\tau)^{3/2}$, we obtain the proportionality,

Here, the ellipses $\dots\sim\int_{0}^{M}ds(\dots)$ 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. for details we note that they behave as $\sim 1/Ms^{1/2}$, matching the above $s^{-3/2}$ law at $s\sim 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 $\tau^{-3/2}$ scaling in this relation is a hallmark of Liouville quantum mechanics. As mentioned in the Introduction, cf. Eq. (14), 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(\beta)\sim\int dk|\langle\phi_{0}|k\rangle|^{2}\exp(-\beta k^{2}/2M)\sim(M/\beta)^{3/2}$, substituting the time arguments, Eq. (74), and continuing to real times, $\tau\to it$, we obtain

Incidentally, the $t^{-6}$ power law scaling coincides with that found for non-interacting spinless fermions .

The application of the same spectral decomposition procedure to the Greens function, $G(\tau)$, yields the result

Using this formula to normalize the OTO function, $F(t)\to F(t)/[G(\beta/2)]^{2}$, we obtain the result describing regime 4 in Eq. (10).

The results derived above can readily be extended to the regime of long times but high temperatures, $\beta\ll M\ll t$. To this end we substitute into Eq. (77) the high temperature expression for the partition sum (12), and use the scaling $G(\beta/2)\sim(J\beta)^{-1/2}$ of the 2-point Greens function to normalize. As a result, Eq. (78) changes to

which now describes the OTO function in the regime 3 of Eq. (10).

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 de-correlations 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\sim N\ln(N)J^{-1}$, where $N$ is the number of particles in the system or, equivalently, a logarithmic measure for the dimension $\sim 2^{N/2}$ of the many-body Hilbert space. This regime is inaccessible to theoretical approaches taking thermodynamic limits $N\to\infty$ 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 $\sim 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\times 3/2$ power law has its origin in the celebrated $t^{-3/2}$ universality of the Liouville quantum mechanics . It is a robust feature of this system that temporal correlations of local operators decay as $\langle\hat{O}(0)\hat{O}(t)\rangle\sim t^{-3/2}$. This power-law may also be interpreted in terms of the $\sqrt{E}$ scaling of the average low-energy many-body DoS, Eq. (13), 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 (47), (48) to be separated by long intervals $\sim|t\pm i\beta/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 ($\gg M$) time intervals involved.

Note added: Recent preprint discusses one-loop (semiclassical) exactness of the partition sum of the Schwarzian theory (12) using the Duistermaat-Heckman formula . It remains to be seen if this technique is also helpful in evaluation of correlations functions.

Acknowledgments: We thank J. Verbaarschot and A.M. García-García for discussions and for sharing unpublished numerical data. A.K. was supported by DOE contract DEFG02- 08ER46482, and acknowledges support by QM2–Quantum Matter and Materials of Cologne University. Work supported by CRC 183 of the Deutsche Forschungsgemeinschaft (project A03).

Appendix A Derivation of Eq. (16)

where $\chi^{a}_{i}$ are $Nn$ independent Grassmann variables, $a=1,\dots,n$ is a replica index, and the replication is introduced to get rid of the partition sum after observables have been computed, $\lim_{n\to 0}Z^{n}=1$. Averaging over disorder, we obtain

where the effective action is given by Eq. (16).

Appendix B Remarks on the Lagrange multiplier

The job of the Lagrange multiplier is to establish, in the best possible way, the constraint (29). It can be understood in analogy $\gamma\leftrightarrow\mu$ to the chemical potential fixing the particle number $\int e^{\phi}\leftrightarrow N$ to a premeditated value $t_{H}\leftrightarrow N_{0}$ in statistical mechanics. Furthering this analogy, fluctuations away from the intended value $\Delta t_{H}[\phi]\equiv\int e^{\phi}-t_{H}$ are obtained as expectation values $\langle\Delta t_{H}[\phi]\rangle_{\phi}=-\partial_{\gamma}\ln Z(\gamma)$, where $Z(\gamma)\leftrightarrow Z(\mu)$ is the functional partition sum depending on the chosen value of $\gamma$, and fluctuations as

Appendix C Gaussian fluctuations

The straightforward Gaussian expansion of the action (33) in $\bar{\varphi}+\delta\varphi$ leads to

The resulting fluctuation determinant is obtained by standard recipes : consider any two independent solutions $z_{1,2}$ of the differential equation

The fluctuation determinant is then given by

where $W=z_{1}^{\prime}z_{2}-z_{1}z_{2}^{\prime}=\pi(1-\lambda)$ is the time-independent Wronskian and $M$ features as the effective inverse ‘Planck-constant’ . Multiplying the exponentiated saddle point action with the factor $D^{-1/2}$ and re-introducing $\beta$ we obtain

which is Eq. (12) up to the normalization prefactor.

Appendix D Computation of the cross correlation Eq. (59)

The central object in Eq. (59) is the expectation value of the fluctuation field, which assumes the form

where $z_{1,2}$ are independent solutions of the differential Eq. (86) satisfying the boundary conditions $z_{1}(-\pi/2)=0$, $z_{2}(\pi/2)=0$ and $W$ is their Wronskian. For times $\tau$ in the bulk of the interval, the boundary regulators $\lambda$ can be neglected and the functions $z$’s entering the fluctuation correlator become

with $W=\pi$. These expressions need to be substituted into Eq. (59) and integrated as

where $\bar{\varphi}=\ln(\cos^{-2}\tau)$ was used. Although this is straightforward in principle, the evaluation is cumbersome in practice. Calculation becomes simpler if we note that (cf. Eq. 48) we may set $\tau_{3}=\tau_{4}+\pi/2$ and $\tau_{1}=\tau_{2}+\pi/2$. Prior to analytic continuation the time arguments lie on the real axis and are ordered as $\tau_{4}<\tau_{2}<\tau_{3}<\tau_{1}$ according to the values of the readout times (48). This structure suggests to define integrals over simplified ‘quadratic’ integration domains,

and apply elementary geometric reasoning to show that Eq. (92) reduces to

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 (60).

After the scaling, $\gamma\to 1$, $\alpha\to\alpha/\gamma$ the Liouville eigenfunctions in the $\phi$-representation are given by Eq. (43) with the replacement $\gamma\to 1+\alpha$. As a consequence, these matrix elements of the operator $\exp(\phi/4)$ assume the form

For generic values of $k,k^{\prime}$ 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

where $K$ is the complete elliptic integral of the first kind. For small values of $\alpha$ the functions $K$ are regular, and for large arguments weakly decay as $\alpha^{-1/4}$. This translates to a logarithmic singularity of the $\alpha$-double integral. What comes to rescue is the convergence generating factor (cf. discussion below Eq. (50)). In the scaled framework, $\gamma\to 1$, this factor is upgraded as $\delta\sim J^{-1}$ to $\delta\gamma\sim 1$ and to the integral we need to perform reads

Appendix F Stationary phase at intermediate times

At intermediate time scales the quench potentials can no longer be neglected and we need to expand around the extrema of the action $S[\varphi]+S_{\alpha_{1}}[\varphi]+S_{\alpha_{2}}[\varphi]$, where $S[\varphi]$ is given by Eq. (33) and $S_{\alpha_{i}}[\varphi]$ by Eqs. (53) with the replacement $\phi\to\varphi$. Presently, we find it convenient to scale the auxiliary variables as $\alpha_{i}\to 2M\alpha_{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 $\bar{\varphi}^{\prime\prime}-2\gamma_{i}\exp(\bar{\varphi})=0$, $\gamma_{0}=1,\gamma_{1}=1+\alpha_{1},\gamma_{2}=1+\alpha_{1}+\alpha_{2},\gamma_{3}=1+\alpha_{2},\gamma_{4}=1$ for the five temporal regimes defined by the time arguments $-\beta/2,\tau_{4},\tau_{2},\tau_{3},\tau_{1},\beta/2$ assume the form $\bar{\varphi}=\ln(\omega^{2}\gamma_{i}^{-1}\cos^{-2}(\omega\tau+\delta))$, where $\omega$ and $\delta$ are free parameters. Requiring the usual boundary singularity $\bar{\phi}(\pm\pi/2)\to+\infty$ and continuity of the solution in the bulk a straightforward fixation procedure yields for these parameters

where Eq. (48) was used. Similarly, the substitution of the stationary solution into the action yields $S[\bar{\varphi}]=-2M\pi\omega_{0}^{2}$. (To derive this result, the regularization of the solution near the boundary by the parameter $\lambda$ of Eq. (35) needs to be taken into account. As before it leads to a cancellation of the formally divergent factor $\gamma 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 (12) and may be ignored. Combining terms, we arrive at Eq. (66).

References