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

TL;DR

This paper analyzes the SYK model's finite temperature properties, revealing a universal power-law decay in out-of-time-order correlations at large times, derived through a Liouville quantum mechanics approach.

Contribution

It introduces an exact integration method over conformal modes to study the SYK model, uncovering a universal power-law decay in correlation functions at large times.

Findings

Out-of-time-order correlation functions exhibit a crossover from exponential decay to a $t^{-6}$ power-law.

The analysis employs a mapping to Liouville quantum mechanics for exact results.

The power-law behavior appears at time scales proportional to the number of particles.

Abstract

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.

Equations

$$\hat{H}=\sum _{ijkl}^N{J}_{ijkl}{\chi }_i{\chi }_j{\chi }_k{\chi }_l,$$

$$M\equiv \frac{N\ln (N)}{64\sqrt{\pi }J}.$$

$$F(t)\equiv \frac{1}{Z{N}^2}\sum _{ij}^N⟨\mathrm{tr}\left(e^{-\frac{\beta }{4}\hat{H}}{\chi }_i(0)e^{-\frac{\beta }{4}\hat{H}}{\chi }_j(t)e^{-\frac{\beta }{4}\hat{H}}{\chi }_i(0)e^{-\frac{\beta }{4}\hat{H}}{\chi }_j(t)\right)⟩,$$

$$t_E\equiv \frac{\ln (MT)}{2\pi T}$$

$$T\gg M^{-1}:$$

$$T\ll M^{-1}:$$

$$Z(\beta )\sim {\left(M/\beta \right)}^{3/2}\exp (2{\pi }^{2}M/\beta ),\beta =\frac{1}{T}$$

$$\rho (E)\propto \theta (E)\sinh \left[2\pi \sqrt{2ME}\right],$$

$$F(t)=\frac{1}{Z{N}^2}\sum _{ij,m_i}⟨⟨m_1|{\chi }_i|m_2⟩⟨m_2|{\chi }_j|m_3⟩⟨m_3|{\chi }_i|m_4⟩⟨m_4|{\chi }_j|m_1⟩e^{-\left(\frac{\beta }{4}+it\right){ϵ}_{m_1}-\left(\frac{\beta }{4}-it\right){ϵ}_{m_2}-\left(\frac{\beta }{4}+it\right){ϵ}_{m_3}-\left(\frac{\beta }{4}-it\right){ϵ}_{m_4}}⟩.$$

$$F(t)\sim |Z(\beta /4+it){|}^4,$$

$$-S[\Sigma,G]={N\over 2}\mathrm{Tr}\log(\partial_{\tau}\delta^{ab}+\Sigma^{ab}_{\tau,\tau^{\prime}})+{N\over 2}\int\!\!\!\int\limits_{-\beta/2}^{\beta/2}d\tau d\tau^{\prime}\left[\frac{J^{2}}{4}\big{(}G^{ab}_{\tau,\tau^{\prime}}\big{)}^{4}+\Sigma^{ba}_{\tau^{\prime},\tau}G^{ab}_{\tau,\tau^{\prime}}\right].$$

$$\displaystyle-(\partial_{\tau}+\Sigma)\cdot G=1;\quad\quad\Sigma=J^{2}\big{[}G\big{]}^{3}.$$

$$G_{\infty }^{ab}(\tau -{\tau }^{\prime })=-\frac{b}{J^{1/2}}\frac{{\delta }^{ab}\mathrm{sgn}(\tau -{\tau }^{\prime })}{|\tau -{\tau }^{\prime }{|}^{1/2}};{\Sigma }_{\infty }^{ab}(\tau -{\tau }^{\prime })=-b^3{J}^{1/2}\frac{{\delta }^{ab}\mathrm{sgn}(\tau -{\tau }^{\prime })}{|\tau -{\tau }^{\prime }{|}^{3/2}},$$

$$G_{\infty }([f],\tau ,{\tau }^{\prime })$$

$${\Sigma }_{\infty }([f],\tau ,{\tau }^{\prime })$$

$$\tau\rightarrow g(\tau):=\tan\bigl{(}\pi\tau/\beta\bigr{)},$$

$$G_{\beta }(\tau -{\tau }^{\prime })=G_{\infty }([g],\tau -{\tau }^{\prime })=-\frac{b}{J^{1/2}}\mathrm{sgn}(\tau -{\tau }^{\prime }){\left[\frac{\pi }{\beta \sin (\pi |\tau -{\tau }^{\prime }|/\beta )}\right]}^{1/2}.$$

$$f(\tau +\beta )-f(\tau )=\beta ;f^{\prime }(\tau +\beta )=f^{\prime }(\tau )\ge 0.$$

$$G_{\beta}([f],\tau,\tau^{\prime})=h^{\prime}(\tau)^{1/4}G_{\infty}\bigl{(}h(\tau)-h(\tau^{\prime})\bigr{)}h^{\prime}(\tau^{\prime})^{1/4},$$

$$⟨{\mathcal{O}}_{\beta }⟩\propto \int \mu [f]\mathcal{D}f{\mathcal{O}}_{\beta }[f]e^{-S_{\beta }[f]}\stackrel{f=g^{-1}\circ h}{=}\int \mu [h]\mathcal{D}h\mathcal{O}[h]e^{-S[h]},$$

$$S[h]=-M\underset{-\beta /2}{\overset{\beta /2}{\int }}d\tau \{h(\tau ),\tau \},$$

$$M(E)=\frac{N\ln (J/E)}{64\sqrt{\pi }J},E>1/M(E).$$

$$h^{\prime }(\tau )\equiv e^{\varphi (\tau )},\varphi (\tau +\beta )=\varphi (\tau ),$$

$$\underset{-\beta /2}{\overset{\beta /2}{\int }}{e}^{\varphi (\tau )}d\tau =t_H,$$

$$S[h]\to S[\varphi ]=\frac{M}{2}\underset{-\beta /2}{\overset{\beta /2}{\int }}d\tau \left[({\varphi }^{\prime }{)}^2-2{\varphi }^{\prime \prime }\right].$$

$$⟨{\mathcal{O}}_{\beta }⟩\propto {Z_{\beta }}^{-1}\underset{\varphi (-\beta /2)={\varphi }_0}{\overset{\varphi (\beta /2)={\varphi }_0}{\int }}\mathcal{D}\varphi \mathcal{O}[\varphi ]e^{-S[\varphi ]},S[\varphi ]=\underset{-\beta /2}{\overset{\beta /2}{\int }}d\tau \left\{\frac{M}{2}\left[({\varphi }^{\prime }{)}^2-2{\varphi }^{\prime \prime }\right]+\gamma e^{\varphi (\tau )}\right\}-\gamma t_H,$$

$$\varphi \equiv \ln \left(\frac{2M}{\gamma }\right)+\phi ,$$

$$S[\phi ]=M\underset{-\pi /2}{\overset{\pi /2}{\int }}d\tau \left(\frac{1}{2}[{\phi }^{\prime }(\tau ){]}^2+2e^{\phi (\tau )}\right),$$

$$\underset{-\pi /2}{\overset{\pi /2}{\int }}{e}^{\phi (\tau )}d\tau =\frac{t_H\gamma }{2M},$$

$$\displaystyle\bar{\varphi}(\tau)\equiv\ln\frac{(1-\lambda)^{2}}{\cos^{2}\Big{(}\tau(1-\lambda)\Big{)}},\qquad\lambda\ll 1.$$

$$\lambda =\frac{8M}{\pi \gamma t_H}.$$

Full text

††journal: Nuclear Physics B

1 Introduction

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

[TABLE]

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 [2, 3, 4, 5] 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) [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 mean-field level the infinite dimensional conformal symmetry gets broken by the interaction self-energy down to the conformal group $\mathrm{SL}(2,R)$ of rational transformations, $t\to\tfrac{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 $\partial_{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^{\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 [22]. 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

[TABLE]

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 [1, 2, 4], $G(\epsilon)\sim|\epsilon|^{-1/2}$ has a divergent amplitude at low energies, the inclusion of Goldstone mode fluctuations shows [21] 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. [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-model — four 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], $\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

[TABLE]

where $Z=\left\langle\mathrm{tr}\left(\,e^{-\beta\hat{H}}\right)\right\rangle$ 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\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 [6, 24] two different regimes characterized by the exponential loss of correlations, depending on whether $t<t_{E}$, or $t>t_{E}$, where

[TABLE]

plays the role of an Ehrenfest time [23] 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.)

[TABLE]

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

[TABLE]

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

[TABLE]

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}$, [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\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’

[TABLE]

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 [20, 13] 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

[TABLE]

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 [21]. 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.

2 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 A)

[TABLE]

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 [3] assume the form of the self-consistent Dyson equation

[TABLE]

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

[TABLE]

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

[TABLE]

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, $\tau\to\tfrac{a\tau+b}{c\tau+d},ad-bc=1$, i.e. the Goldstone mode manifold is identified as the group of time-reparameterizations modulo the unbroken $\mathrm{SL}(2,R)$-transformations.

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]$ [27]: the reparameterization

[TABLE]

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

[TABLE]

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

[TABLE]

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

3 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. [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:[-\beta/2,\beta/2]\to[-\beta/2,\beta/2]$ mapping the finite imaginary time interval onto itself (cf. Eq. (23)), we find it convenient to work with functions $h:[-\beta/2,\beta/2]\to\mathbb{R}$ which have the full real numbers as their image. (Within the $[-\beta/2,\beta/2]\to[-\beta/2,\beta/2]$ framework, the weak symmetry breaking part of the effective action picks up potential contributions under re-transformations which are unpleasant to work with.) The passage to this description is defined by $h(\tau)\equiv g(f(\tau))$, where $g$ is the $\tan$-function, Eq. (21). The function $h$ then necessarily has a singularity at some $\tau^{*}\in[-\beta/2,\beta/2]$, where $f(\tau^{*})=\pm\beta/2$ (cf. Fig. 2.) A shift of the time axis, which is part of the $\mathrm{SL}(2,R)$ invariance manifold, can be applied to shift the singularity to the boundaries of the interval $\tau^{*}=\pm\beta/2$.

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

[TABLE]

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,

[TABLE]

where we omitted the subscript ${}^{\prime}\infty^{\prime}$ for brevity. Here, $\mu[h]\propto{\prod_{\tau}}[h^{\prime}(\tau)]^{-1}$ [21] is the group measure whose reparameterization invariance implies $\mu[f]\mathcal{D}f=\mu[h]\mathcal{D}h$. The action $S[h]$ in Eq. (25) describes the cost of $h$-fluctuations due to the presence of the derivative $\partial_{\tau}$ in Eq. (16). Requiring invariance under the $\mathrm{SL}(2,R)$-reparameterization group and assuming locality in time, one may argue [4] that this action must be the integrated Schwarzian derivative

[TABLE]

where $\{h,\tau\}\equiv(\tfrac{h^{\prime\prime}}{h^{\prime}})^{\prime}-\tfrac{1}{2}(\tfrac{h"}{h^{\prime}})^{2}$ and $M$ is a coupling constant of dimensionality $[\mathrm{time}]$111 The action $S_{\beta}[f]$ is obtained from here as $S_{\beta}[f]\equiv S[g\circ f]=S[f]-\frac{2\pi^{2}M}{\beta^{2}}\int\limits_{-\beta/2}^{\beta/2}d\tau\,f^{\prime 2}$, cf. Ref [4]. . On dimensional grounds, this constant determines the time scale, $t_{*}\sim 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 (16) 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

[TABLE]

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

[TABLE]

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$ [21]. 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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

4 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

[TABLE]

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

[TABLE]

the constraint becomes

[TABLE]

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 [27] 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,

[TABLE]

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

[TABLE]

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

Now is a good time to discuss the status of the Lagrange multiplier, $\gamma$. This parameter was introduced to enforce the constraint, Eq. (34). In conventional variational calculus, the value of a Lagrange multiplier is fixed by requiring that the constraint is satisfied on the variational configuration $(\bar{\varphi})$. Presently, however, we have one more free variable in the picture, viz. the diverging boundary values $\varphi_{0}$. We may use this freedom to make sure that for any given value of $\gamma$ the constraint is satisfied (we saw above that this is achieved by the value $\exp(\varphi_{0})=(2/\pi\lambda)^{2}=(\gamma t_{H}/4M)^{2}$.) This gives us the freedom to fix $\gamma$ at will, and the best choice is $\gamma=J$, i.e. the largest possible value of a variable with dimension $[\mathrm{energy}]$ in the problem. A large value of the Lagrange multiplier is favored because it minimizes fluctuation-violations of the constraint (much like particle number fluctuations in statistical mechanics are minimized by large chemical potentials. For a few more remarks substantiating this point, see B.) Our rationale, thus, is to impose a fixation $\gamma=J$ throughout. This determines the boundary values of the Liouville quantum mechanics which in the language of the original variable, $\phi$ read

[TABLE]

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

[TABLE]

independent of the regularization.

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.

5 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. (31) lends itself to a straightforward quantum mechanical interpretation. To this end, consider the Liouville Hamiltonian operator

[TABLE]

and the corresponding imaginary time propagator:

[TABLE]

Comparison with Eq. (31) shows that the regularized Schwarzian partition function can be identified with the quantum propagator

[TABLE]

provided we have the identification

[TABLE]

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

[TABLE]

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

[TABLE]

Into this expression we substitute $\phi_{1}=\phi_{2}=\phi_{0}$, which means that we have to consider the (real) wave function amplitudes

[TABLE]

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

[TABLE]

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. [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 $\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.

6 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

[TABLE]

where the angular brackets denote both, disorder averaging, and the quantum thermal average $\mathrm{tr}(\exp(-\beta\hat{H})(\dots))$, and $T_{\tau}$ 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)

[TABLE]

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 $\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,

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

and $S[\phi]$ is given in Eq. (40).

Eq. (52) and (33) define the path integral we need to compute. Its exponent suggests an interpretation in terms of an effective quench dynamics:222 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. [30] 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.

7 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

[TABLE]

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,

[TABLE]

where we defined

[TABLE]

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

[TABLE]

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$,

[TABLE]

In addition, the fluctuation offset $\delta\varphi$ appears in the pre-exponential factors, and in the quench potentials $\alpha_{i}e^{\varphi}$. 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(\delta\varphi(\tau_{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 $S_{\alpha_{i}}[\bar{\varphi}+\delta\varphi]=S_{\alpha_{i}}[\bar{\varphi}]+\sum_{n=1,2}S_{\alpha_{i}}^{(n)}[\bar{\varphi},\delta\varphi]$, where terms with superscript $(n)$ are of $n$th order in $\delta\varphi$, produces two contributions, $\langle S_{\alpha_{i}}^{(2)}\rangle-\tfrac{1}{2}\langle(S_{\alpha_{i}}^{(1)})^{2}\rangle$, where the angular brackets denote functional averaging over the fluctuation action (58). These lead to a weak renormalization of the $i$th Green function $G_{\beta}$ 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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

8 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}$:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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. [6] which predicted a decay $\sim t\exp(-\pi t/\beta)$.

For completeness we mention that the above discussion assumed a fixation of the real parts $\mathrm{Re}(\tau_{i})$ of the time arguments relative to that, $\tau^{\ast}=\beta/2$, of the singularity which is necessarily present in a bijective transformation $[-\beta/2,\beta/2]\to\mathbb{R}$ (cf. discussion in the beginning of section 3.) This fixation may be abandoned by introducing a shift parameter $\tau_{i}\to\tau_{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 $\tau_{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:

9 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$.

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 (48) at $it=\tau\in\mathbb{R}$. Here, $\tau$ serves as a real time parameter which will eventually be continued into the complex plane, and $s\in[0,\beta/4]$ is the parameter controlling the relative position of the observation times and the time, $\tau^{\ast}$, of the reparameterization singularity (cf. discussion at the end of the previous section). The correlation function is then obtained as

[TABLE]

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

[TABLE]

where $\tau_{ij}=\tau_{i}-\tau_{j}$ are the time differences

[TABLE]

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

[TABLE]

Here, we used that for the diverging boundary conditions, $\phi_{0}\to\infty$, the boundary wave functions, $\langle\phi_{0}|k\rangle=\mathcal{N}(k)K_{2ik}(2\sqrt{2M}e^{\phi_{0}})\to\mathcal{N}(k)\times\mathrm{const.}$, 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 $|\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,

[TABLE]

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. [21] 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

[TABLE]

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(\tau)$, yields the result

[TABLE]

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

[TABLE]

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

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 [26]. 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 [32] discusses one-loop (semiclassical) exactness of the partition sum of the Schwarzian theory (12) using the Duistermaat-Heckman formula [33]. 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)

In this appendix, we briefly review how the action (16) emerges as an effective description of the system described by the Hamiltonian (1). To this end, consider a Grassmann coherent state functional integral representation of the replicated partition sum $Z^{n}\equiv(\mathrm{tr}(\exp(-\beta\hat{H}))^{n}$,

[TABLE]

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

[TABLE]

where $\tilde{G}^{ab}_{\tau,\tau^{\prime}}=-N^{-1}\sum\limits_{i}^{N}\chi_{i}^{a}(\tau)\chi_{i}^{b}(\tau^{\prime})$. To pass from the abbreviation, $\tilde{G}$, for a Grassmann bilinear to an integration variable, $G$, we introduce

[TABLE]

in the functional integral, which entitles us to replace $\tilde{G}\to 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 $\left|\partial_{\tau}\delta^{ab}+\Sigma^{ab}_{\tau,\tau^{\prime}}\right|^{N\over 2}$. We re-exponentiate the latter to obtain the representation

[TABLE]

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

[TABLE]

The evaluation of these formulas on the partition sum (89) yields $\langle\Delta t_{H}[\phi]\rangle_{\phi}=\gamma^{-1}$ and $\mathrm{var}(\Delta t_{H}[\phi])=\gamma^{-2}$, respectively. This shows that large values of $\gamma$ increase the accuracy of the formalism. Our logics is to implement the constraint through freely adjustable parameters (notably the values $\phi_{0}$ of the field at the boundaries of the temporal domain) and choose for $\gamma$ the largest value consistent with the low energy nature of the theory, $\gamma\sim J$. For these values, fluctuations in the constraint variables, which have dimensionality $[\mathrm{time}]$ are of $\mathcal{O}(J^{-1})$ which is the smallest time scale resolvable in the theory.

Appendix C Gaussian fluctuations

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

[TABLE]

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

[TABLE]

For later convenience, we choose

[TABLE]

The fluctuation determinant is then given by [34]

[TABLE]

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’ [34]. Multiplying the exponentiated saddle point action with the factor $D^{-1/2}$ and re-introducing $\beta$ we obtain

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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,

[TABLE]

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

[TABLE]

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

Appendix E Computation of the matrix elements $\langle k^{(\alpha)}|e^{\frac{\phi}{4}}|k^{{}^{\prime}(\alpha^{\prime})}\rangle$

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

[TABLE]

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

[TABLE]

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

[TABLE]

where $C$ is a constant of order unity.

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

[TABLE]

where $\omega$ is specified in Eq. (66). (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, $|\mathrm{Im}(\tau_{i})|=t/2$, and an approximation $\varphi=\ln(\omega^{2}\gamma_{i}^{-1})\pm 2i(\omega\tau-\delta)$ is possible, where the sign is postive/negative for $\tau_{4,3}$/$\tau_{1,2}$.) The substitution of these configurations into the pre-exponential factor yields

[TABLE]

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

**