On the replica structure of Sachdev-Ye-Kitaev model

We investigate existence of replica off-diagonal solutions in the field-theoretical description of Sachdev-Ye-Kitaev model. To this end we evaluate a set of local and non-local dynamic correlation functions in the long time limit. We argue that the structure of the soft-mode Schwarzian action is qualitatively different in replica-diagonal vs. replica-off-diagonal scenarios, leading to distinct long-time predictions for the correlation functions. We then evaluate the corresponding correlation functions numerically and compare the simulations with analytical predictions of replica-diagonal and replica-off-diagonal calculations. We conclude that all our numerical results are in a quantitative agreement with the theory based on the replica-diagonal saddle point plus Schwarzian and massive Gaussian fluctuations (the latter do contain replica off-diagonal components). This seems to exclude any contributions from replica-off-diagonal saddle points, at least on the time scales shorter than the inverse many-body level spacing.

The model is represented by interacting Majorana fermions with quenched random matrix elements. As such, it naturally admits a description in terms of a replica field theory [52]. The structure of this theory in the replica space has important bearings on all aspects and applications of the SYK-like models. In particular, an existence and JHEP09(2019)057 properties of a glass phase is most naturally discussed in terms of the replica symmetry breaking (RSB). Indeed, RSB was first introduced by Parisi [53,54] to describe glass transition in the Sherrington-Kirkpatrick model [55]. The existence of the glass phase in SYK and similar models is a subject of intense discussions since the very introduction of the model [1,44,[56][57][58]. A recent discussion, ref. [59], came on the side of the absence of the glass phase.
Non-trivial replica structures, associated with some form of RSB, were discussed in many other fields of the statistical physics. Most relevant to the present context are replica studies of the random matrix theory. The latter may be classified as SYK 2 model (as opposed to SYK 4 , discussed in this paper). The corresponding replica field theory of SYK 2 is known as non-linear sigma model. Its long time (i.e. small energy) correlation functions were understood in terms of the broken replica symmetry [60][61][62][63][64]. To some extent these studies mirror Altshuler-Andreev description [65] in terms of the broken supersymmetry. The important point is that RSB is only noticeable on the time scale associated with the inverse level spacing (Heisenberg time) and practically does not have consequences at shorter times.
It was recently suggested [66] that non-trivial replica structure (i.e. replica non-diagonal saddle point) may be responsible for the behavior of the structure factor of SYK 4 model at a time scale parametrically shorter than the Heisenberg time. Another recent study [67] discusses thermodynamic relevance of the replica off-diagonal saddle points in SYK model.
The goal of this paper is to investigate possible signatures of RSB and replica offdiagonal saddles on the behavior of correlation functions at moderately long times. By those we understand time scales longer than τ > N/J, yet shorter than the Heisenberg time (i.e. inverse many-body level spacing). The correlation functions considered here are motivated by mesoscopic physics [68], where one is interested in quench disorder averages of higher moments of certain quantum observables. (One may also look for an entire probability distribution function of a given observable over quench disorder realizations.) Here we show that the corresponding correlation functions exhibit qualitatively distinct behavior being calculated on replica diagonal vs. replica off-diagonal saddle points. The difference stems from the ways the corresponding saddle points break the reparameterization symmetry [2,3,69] of the model. The distinct patterns of the symmetry breaking are reflected in the structure of the low-energy (Schwarzian) action. We found that in case of replicadiagonal saddle points the latter consists of n (number of replica) independent Schwarzians. However, for a generic replica off-diagonal saddle point there is only one Schwarzian degree of freedom, while the remaining n − 1 acquire a stiffer action. These observations translate into a different behavior of correlation functions at moderately long times.
We then perform a detailed comparison of our analytical expectations with numerical simulations of N ≤ 32 SYK 4 model. Our simulations use exact diagonalization and exact matrix elements to evaluate corresponding "mesoscopic" correlation functions. We consider p-th moments, p = 1, . . . , 5, of both site-local and site-non-local two-point correlation functions. The comparison shows no evidence for contributions from replica off-diagonal saddles. On the contrary, all the data may be quantitatively accounted for by the theory based on replica diagonal saddle point along with reparameterization fluctuations and JHEP09(2019)057 massive Gaussian fluctuations around it. The massive fluctuations, which include replica off-diagonal components, must be retained to account for small site non-local correlations.
This work is organized as follows: in section 2 we review SYK model and its description in terms of the replica field theory. In section 3 we discuss soft reparameterization modes and how their action is different between replica diagonal and replica off-diagonal saddle point configurations. In section 4 we discuss consequences of these differences for the longtime behavior of mesoscopic correlation functions. We put these differences to numerical test in section 5. In section 6 the similar program is implemented to a different family of mesoscopic correlation functions -those with site non-local correlations. We present brief conclusions in section 7. A number of technical details are relegated to appendices A-F.

Model and definitions
The real-fermion version of SYK 4 model is formulated in terms of Majorana fermions χ i on lattice sites i. It is determined by the Hamiltonian where matrix elements J ijkl are random statistically independent Gaussian distributed variables with zero mean and a variance given by The standard field theoretical treatment of the SYK model [46,69] employs the replicatrick. It replicates the fermionic degrees of freedom as χ i → χ a i , where a = 1, . . . n, allowing for a direct averaging over the random couplings J ijkl . To arrive at an effective bosonic field theory, describing the behavior of the model at low energies and long times, one then integrates the fermionic degrees of freedom, by introducing a replica-matrix valued field Eq. (2.3) is enforced by inserting the functional δ-function in the replicated partition func- where the matrix field Σ ab τ τ is the Lagrangian multiplier dual to G ab τ τ . After integration of the fermionic degrees of freedom, one arrives at the following action In the large N limit, where N is the number of sites, the properties of the model are determined by the saddle point of the path integral over the effective fields G ab τ τ and Σ ab τ τ . The corresponding saddle point equations read . At this junction the standard choice [6,[45][46][47]49] is to look for a replica-diagonal saddle point solution of the formĜ = δ ab G τ τ and correspondinglŷ Σ = δ ab Σ τ τ . Let us emphasize that, although we call such a choice replica-diagonal, fluctuations around the replica-diagonal saddle may and should include replica-off-diagonal components δG ab τ τ . We discuss them in detail in appendix B. In the long time limit one can neglect the time-derivative term in the second equation in (2.6), which makes the equations scale invariant. In such a limit one may look for replica-off-diagonal solutions in a separable form, i.e. where matrix form in replica and time spaces separates as:Ĝ = g ab G τ τ ,Σ = σ ab Σ τ τ . (2.7) Here G τ τ and Σ τ τ are traditional replica-diagonal solutions, while time-independent symmetric n × n matrices g and σ satisfy: This particular form is motivated by the fact that it allows to naturally keep the conformal structure of the long-time effective theory [2,3] -the feature that was proven to be central to the physics of the SYK model. Combining the ansatz (2.7) with the usual replica-diagonal conformal solution, one finds: where replica matrices g and σ satisfy eqs. (2.8).
The n × n, where n is number of replica, matrix equations (2.8) admit a wealth of both diagonal and off-diagonal solutions. It is thus necessary to spell out selection criteria on which of these solutions should be taken into account and why. The most natural of such criteria seems to be a requirement of having a minimal action (i.e. free energy). In particular one may ask if the widely accepted choice g = σ = δ ab indeed has the smallest action. In appendix A we show that one can find a discrete set of solutions of the form g ab =gδ ab + g(1 − δ ab ), (2.10) whereg and g are n-dependent complex numbers. Moreover, in the n → 0 limit the (real part of) corresponding action is smaller than that on the diagonal (i.e.g = 1, g = 0) solution. Similar conclusions were recently reached in ref. [67]. The question thus arises whether these (or others) replica-off-diagonal solutions are indeed relevant for the physics of the model. This question is farther complicated by the fact that besides the saddle point action one needs to evaluate fluctuation determinants and perform summation over the set replicaoff-diagonal saddles for any desired observable. Since we do not know how to perform this program in general, we seek for generic signatures, which help to distinguish between diagonal and off-diagonal solutions. Below we argue that long time behavior of certain correlation functions serves as a sensitive test for the presence of the off-diagonal components.

JHEP09(2019)057
To argue why this is indeed the case one needs to consider a structure of soft-mode fluctuations around diagonal and off-diagonal solutions. For the conformal solutions of the form eq. (2.9) such soft modes are given by reparameterization fluctuations [2,3,6,[45][46][47]49].

Reparameterization fluctuations
In the conformal limit (i.e. neglecting ∂ τ term) the action (2.5) and the saddle point equations eqs. (2.6), are invariant under the time reparametrization transformations where G ab (τ 1 , τ 2 ) and Σ ab (τ 1 , τ 2 ) are conformal solutions (2.9). Here f a (τ ) with a = 1, . . . n is a replica-specific reparametrization transformation. This defines the symmetry group G of the action (2.5) in the infra-red limit, G = ⊗ n a=1 Diff(R), where Diff(R) denotes the diffeomorphism group of time axis. The product over replicas reflects the fact that reparametrization transformations can be chosen independently in different replicas, i.e. f a (τ ) = f b (τ ). The symmetry under time-reparametrizations is a crucial property, that relates the SYK model to the AdS 2 gravity theories [4][5][6][7][8][9][10][11][12]70]. This time-reparametrization symmetry is however spontaneously broken by the saddle point solutions (2.9) down to the subgroup H = SL(2, R) ⊂ G resulting in the appearance of a soft modes, which span the coset G/H. Specifically, the group H is formed by all Möbius maps of the form h(τ ) = (Aτ + B)/(Cτ + D) with AD − BC = 1.
Here, the major difference shows up between the diagonal and off-diagonal saddle point solutions. For the diagonal case the subgroup is H = ⊗ n a=1 SL(2, R). Indeed, for a diagonal saddle point the independent Möbius maps h a (τ ) = (A a τ + B a )/(C a τ + D a ) may be taken for each replica, leaving the diagonal solution (2.9) invariant. The diagonal soft mode coset is thus G/ H = ⊗ n a=1 Diff(R) / ⊗ n a=1 SL(2, R) . This should be contrasted with the off-diagonal case, where the subgroup is H = SL(2, R) -the same for all replicas. Indeed, performing different Möbius transformations in different replicas does not leave (2.9) invariant, if g and σ have off-diagonal components. 1 Therefore the coset is different: G/H = ⊗ n a=1 Diff(R) /SL(2, R) (similar structure of the coset space was mentioned recently in ref. [71] in context of two-boundary Jackiw-Teitelboim Gravity). The different structure of the coset is reflected in the soft mode action.
The latter action originates from the explicit breaking of the reparametrization symmetry by the time derivative term δ ab ∂ τ . In the diagonal case, where the coset is the product of n independent Diff(R)/SL(2, R) components, the corresponding action is the sum of n Schwarzian derivatives [69] where M ∼ N log N is the mass of the soft fluctuations [46].

JHEP09(2019)057
In the off-diagonal case the subgroup H consists of a single SL(2, R), suggesting that only a single degree of freedom is governed by the Schwarzian action. Indeed, the explicit calculation, outlined in details in appendix E, shows that the off-diagonal matrix elements g ab = 0 in the saddle point solution generate additional terms in the action for reparametrization fluctuations. These terms overpower n − 1 Schwarzian derivatives in the long-time limit.
Let us use the exponential representation of reparametrizations [46] f a (τ ) = τ exp[φ a (τ )]dτ, (3.4) which has an advantage that the corresponding invariant integration measure is flat in φ a variables. In this representation, the additional action can be cast in the form of an effective potential: Here the integration goes along the line C = (τ 1 (τ ), τ 2 (τ )) drawn in R 2 space of two times, at which two reparametrizations take equal values, f a (τ 1 (τ )) = f b (τ 2 (τ )). When expanded in small deviations φ a − φ b , each term in the action eq. (3.5) acquires the form of a "mass" term where in the last expression the integral already goes along the straight line. It is clear that this term is minimized when reparametrizations in all replicas are identical and penalizes deviations from such configuration. To formalize this observation we introduce new variables as φ a = Φ + ϕ a , where n a=1 ϕ a ≡ 0 and therefore Φ = 1 n n a=1 φ a . Then the soft mode action for, e.g., off-diagonal ansatz (2.10) takes the form where we have used that a =b (φ a − φ b ) 2 = ab (ϕ a − ϕ b ) 2 = 2n a ϕ 2 a , since a ϕ a = 0. In the long time limit the last term here suppresses fluctuations of n − 1 degrees of freedom ϕ a , leaving the single degree of freedom Φ, to be governed by the Schwarzian action. This effectively locks reparameterization degrees of freedom in different replicas to Finally, let us mention that the structure of action eq. (3.7) is consistent with the coset space G/H of the replica off-diagonal SYK action. For the infinitesimal reparametrizations f a (τ ) = τ + a (τ ) the phases φ a (τ ) a (τ ). We see that the action eq. (3.5), if written in terms of a (τ ), remains massless vis-a-vis n degrees of freedom. However only single degree of freedom is "super soft": The locking of reparameterization modes in different replicas, eq. (3.8), for off-diagonal saddle points has important consequences for long time behavior of the correlation functions, which we explore in the next section.

Site-local correlation functions
It is well known, that the reparameterization fluctuations modify the long-time decay of correlation functions. The simplest example is the two-point site-local function: While at short times, 1/J < |τ | < N/J, the decay is governed by the conformal mean field behavior G ∼ |τ | −1/2 , eq. (2.9), its long-time behavior, |τ | > N/J, is very different: G ∼ |τ | −3/2 due to the effect of the reparametrization fluctuations [46]. Moreover, the 2p-point correlation functions (p < N ) of the form at long time decay with the same universal exponent −3/2, i.e. G 2p ∼ |τ | −3/2 [46]. The short time behavior is, of course, p-dependent: G 2p ∼ |τ | −p/2 . It is important to notice that the angular brackets in eqs. (4.1) and (4.2) imply both quantum mechanical groundstate expectation value (hereafter we restrict ourselves to zero temperature) along with the averaging over disorder realizations. We now introduce different objects, inspired by mesoscopic fluctuations physics [68] [ where |GS stays for a disorder specific ground-state of the SYK 4 model (the same for all p expectation values), while . . . dis denotes averaging over realizations of random matrix elements J ijkl . In the replica formalism, the correlation function eq. (4.3) can be written as where angular brackets denote averaging with respect to the replicated action and all the replicas a 1 , . . . , a p are different. The leading contribution to the correlation function eq. (4.4) is given by the product of replica-diagonal contractions. Indeed, each contraction of fermions with different replicas enforces the equality of the sites of the contracted fermions, for example χ a 1 i 1 χ a 2 i 2 ∝ g a 1 a 2 δ i 1 i 2 thus eliminating one summation over sites. Such

JHEP09(2019)057
contribution is therefore suppressed by the factor 1/N (in case of the replica diagonal saddle point, g a 1 a 2 = 0 and such contractions originate from Gaussian fluctuations of δG a 1 a 2 and δΣ a 1 a 2 , bringing additional factors of 1/N ). The leading contribution from the product of replica diagonal contractions has furthermore to be averaged over the reparametrization fluctuations where G aa (f a (τ ), f a (0)) is given by eqs. (3.1), (3.4).
In the case of the off-diagonal saddle point, the reparametrizations are locked, eq. (3.8), and therefore the integration in eq. (4.5) runs over the single field Φ. This makes eqs. (4.2) and (4.4) essentially equivalent in the long time regime. One thus expects to find On the other hand, in the replica diagonal case the reparametrizations are not locked, the integration in eq. (4.5) runs over p independent fields and one expects [G(τ, 0)] p dis ∼ τ −3p/2 again in the long time regime. For short times reparameterizations are not relevant and one expects mean-field [G(τ, 0)] p dis ∼ τ −p/2 irrespective of the replica structure. To summarize: This can be checked numerically to distinguish between diagonal and off-diagonal scenario.

Numerical results for site-local correlation functions
The basic quantity for numerical calculations is the two-time ground-state expectation value: In the second equation |n denote many-body excited states (with the parity opposite to that of the ground-state). Numerically, the correlation function eq. (5.1) is calculated from the spectrum of energies and matrix elements obtained by exact diagonalization (see appendix D for details). The correlation function eq. (5.1) is then used to construct the higher order correlation functions as defined by eq. (4.3). Numerical results for the correlation function eq. (4.3) are shown in figure 1. The correlation function G(τ, 0) dis (p = 1) is shown in the left panel in figure 1. Its time-decay exhibits three qualitatively different regimes. At short times (1 τ 10 in units of 1/J) the correlation function decays as τ −1/2 . This behavior corresponds to a saddle point solution, eq. (2.9). At longer times, for 10 τ 100, the time decay changes to τ −3/2 . Such a behavior signals the dominant effect of soft reparametrization fluctuations around the saddle point, as described in ref. [46]. At still longer times, τ 100, the time decay of the correlation function is dominated by a first excited many-body state |n = 1 JHEP09(2019)057 (we'll refer to it as "two-level" system), due to the discreteness of the energy spectrum in a finite size system. The crossover to the two-level regime at long times is quantified on the right panel of figure 1. In that panel, the dashed lines correspond to the calculation of the correlation functions taking into account the two lowest energy levels, |GS and |n = 1 , only.
The right panel in figure 1 shows that the correlation functions ( [G(τ, 0)] p dis ) 1/p , calculated for different p, coincide in a wide time range, which includes both mean-field and reparameterization dominated regimes. Comparing this behavior with the theoretical expectations, eq. (4.6), we conclude that it is consistent only with the replica-diagonal structure of the saddle point. We present an additional independent support to this conclusion by considering site non-local correlation functions in the next section.
Eventually graphs for different p diverge on approaching the two-level system regime. This latter behavior may be quantitatively explained assuming some (independent) distribution functions for matrix elements GS|χ i |1 and energy splitting E 1 − E GS (notice that since the ground-state and the excited state belong to different parity sectors, there is no repulsion between them). See appendix F for more details on the two-level regime. To the best of our knowledge, it is not known how to incorporate two-level regime into the replica filed-theory discussed here (see ref. [16] for an alternative approach). The situation is very different in SYK 2 model, where the corresponding filed-theory is rotationally invariant in the replica space, allowing for the treatment of RSB at the two-level energy scale [60][61][62][63][64].

Site non-local correlation functions
The existence of the replica off-diagonal solutions may be also detected by considering site non-local correlation functions of the type: with i = j. The advantage of this object is that it vanishes, being calculated at the replicadiagonal saddle point (without account for massive fluctuations), but does not vanish, being calculated at the off-diagonal saddle point. To see this we rewrite it in the replica formalism as, where in the second line we disregarded Gaussian fluctuations and used the site-locality of the saddle point correlation functions (both replica diagonal and off-diagonal ones). Since all replica indexes a 1 , . . . , a 2p are distinct here, it is clear that the second line in eq. (6.2) is zero on the diagonal saddle point. To estimate it in the replica non-diagonal saddle point we consider block-diagonal matrices g and σ, consisting of n blocks each of the size 2p × 2p along the main diagonal (see figure 2). In the replica-limit n → 0, one remains with saddle point equations for a single p × p block of the matrix in figure 2. Here we perform explicit calculations for the case p = 1. Using the saddle point ansatz eq. (2.9), we obtain the correlation function D 2 in the form where we use the saddle point matrix g consisting of 2 × 2 blocks. Replica non-diagonal solutions of eqs. (2.8) for 2 × 2 blocks read Eq. (6.3) describes the behavior of the correlation function at short times, when the influence of reparameterization fluctuations is negligible. To obtain the correct time dependence at longer times, replica off-diagonal correlation function has to be averaged over the

JHEP09(2019)057
reparametrization fluctuations (see appendix C for details). Since the reparameterizations are locked according to eq. (3.8), one integrates over a single reparameterization degree of freedom for the both replicas involved in D 2 (τ ). This leads to: We come back now to the replica diagonal scenario, which leads to vanishing result for the non-local functions (6.1), being calculated at the diagonal saddle point. However one can include massive (with the mass of order N ) Gaussian fluctuations δG ab and δΣ ab around the diagonal saddle point to find a non-zero result for the first line in eq. (6.2) (see appendix B for detailed derivation). This leads to: for τ < N/J. Subsequent averaging over the reparameterization fluctuations around the replica-diagonal saddle point with two independent (unlocked) reparameterization modes, one for each replica, results in: Analytical calculations of the correlation function eq. (6.2) for p > 2 result in the following general relation between the correlation functions for different powers p  As one can see from eqs. (6.3), (6.6), for short times the two above mentioned scenarios differ only in the scaling of the correlation function D 2 (τ ) with the number of sites N , while for long times both the scaling with N as well as the predicted time dependence become different. Therefore, the time dependence as well as the dependence on the total number of sites N can be used do discriminate between eqs. (6.5) and (6.6).
Results of numerical calculations of the dependence D 2p (τ ) are shown in figures 3 and 4. Figure 3 shows the time dependence of the correlation function D 2p (τ ) for different p, for our largest system, N = 32. First, one notices the non-monotonous dependence of the correlation functions on time. This short time behavior originates from the fact that equal time expectation χ i (0)χ j (0) = 0 for i = j for N = 8, 16, 24, 32, . . ., which belong to the orthogonal symmetry class [14]. Indeed, from anticommutation of Majoranas one concludes that χ i χ j is pure imaginary. On the other hand, for orthogonal symmetry classes, there is a representation of Majorana operators with all matrix elements n|χ i |m being real.

JHEP09(2019)057
a qualitative agreement with the field theory. However, while p = 1 function is consistent with eqs. (6.6), (6.7), the p ≥ 2 functions exhibit fast crossover to the two-level regime. We thus are not able to verify the time dependence of eqs. (6.8) even for our largest system of N = 32 for p ≥ 2.
We can, however, verify the N -dependence of eqs. (6.8). Numerical results for the dependence of the correlation functions on the number of sites N are shown in figure 4. One can see that the correlation functions rapidly decrease with increasing N while keeping qualitatively the same time dependence. This is in accord with the predictions from the Gaussian fluctuation expansion around the replica-diagonal saddle point. The best fit, see inset in figure 4, for power law dependence on N is N −3.36 , which is close to the power N −3 , following from the Gaussian fluctuation expansion. In contrast, the replica off-diagonal saddle point predicts no suppression of the site off-diagonal correlation functions with N , eq. (6.8). Once again we conclude that the numerics is consistent with the replica-diagonal theory and is inconsistent with the off-diagonal saddle points.

Conclusions
We have examined signatures of the replica off-diagonal saddle points in the field theory treatment of the SYK 4 model. We have argued that such off-diagonal elements affect the coset manifold G/H of reparameterization soft modes and thus change the expected longtime behavior of "mesoscopic" correlation functions. Comparing to numerically evaluated corresponding correlation functions for N ≤ 32 SYK 4 model, we conclude that they do not show any evidence for replica off-diagonal saddle points. On the contrary, all correlation functions (both site local and non-local) are in a good agreement with the expectations stemming from the replica diagonal saddle point plus Gaussian fluctuations. The latter do include replica off-diagonal components, of course.
We conclude thus that we do not detect any evidence for replica off-diagonal saddle points, at least for time scales shorter than the inverse spacing between the ground-state and the first many-body excited state. At longer times the data is well described by statistics of the two-level system. We stress, though, that such two-level description is outside of the field theoretical treatment, we base our conclusions at. It remains to be seen if our conclusions can be reconciled with the results of refs. [66,67]. On the other hand, our findings are in line with no evidence for glassy behavior in the SYK model reported in ref. [59].
Does it mean that there is no room for replica symmetry breaking and replica offdiagonal structure in SYK-like models? In our opinion such conclusion is premature. One may investigate deformed models, such as eg. SYK 4 + SYK 2 [39][40][41]. Our preliminary investigation [72] points to similarities between level and eigenfunction statistics of such models to those of random regular graphs (RRG) [73][74][75] (see also [76,77] for the related studies). In the case of RRG replica symmetry breaking was argued to be a proper framework to describe an observed phenomenology [75]. One is thus justified to expect that phenomenologically similar deformed SYK models may admit similar replica symmetry broken description. However, one should probably conclude that the undeformed In the replica-limit n → 0, equation for z can be written as (z − 1)(z 2 + z + 2) = 0 (A.5) The solution z = 1 meansg = g, which cannot satisfy eq. (A.2). The solutions of equation result in the nontrivial replica non-diagonal saddle-points that are discussed below. For z = −1+i √ 7 2 , we obtain , as well as for the replica diagonal saddle pointg = 1/ √ J, g = 0. Therefore, the difference between the action at the replica-diagonal and at the replica-non-diagonal saddle points comes from the logarithmic terms only. It is given by The solution z = 0 results in the completely replica off-diagonal matrix g, which clearly contradicts numerical results in section 5. From the solution z = i it followsg = ig. Substituting this relation in eq. (A.2) for p = 2, we obtain, as one possible solution, the matrix g as given in eq. (6.4).

B Fluctuation expansion around the replica-diagonal saddle point
To obtain the action for massive fluctuations around the replica diagonal saddle point, we adopt the ansatz Using the symmetry of the fermionic Green functions with respect to exchange of arguments we rewrite the quadratic action in the form Correlation functions of the fields δG and δΣ are obtained by taking derivatives, for instance δΣ ab τ 1 τ 2 δΣ ab τ 3 τ 4 (δG ab τ 5 ,τ 6 ) 4 = 2 6 δ 6 e −δS ab (2) δh ab τ 1 τ 2 δh ab τ 3 τ 4 (δj ab τ 5 τ 6 ) 4 (B.10) Integrating the fields δΣ, δG, we obtain From eq. (B.11) we read off the following nonzero contractions of the fields δG and δΣ δG ab The additive form of the fluctuation action eq. (B.3) implies the multiplicative form of the correlation function eq. (6.2) as a product over the pairs of replicas (B.14)
Here (2p − 1)!! denotes the number of ways to assign pairs of replicas. The fluctuation expansion around the replica-diagonal saddle point results in The same conclusion can be made by means of 1/N diagrammatic expansion without employing the replica-trick (see figure 5) [45].

C Correlation function in the replica non-diagonal saddle point
The replica non-diagonal saddle point assigns a non-zero average to the correlation function χ a i (τ )χ b i (τ ) . Calculation of high-order correlation functions, such as the one defined by eq. (6.1), requires building of contractions of pairs of fermions corresponding to the nonzero average in the saddle point, which leads to the general form as given by eq. (B.14), but with the function D 2 calculated in the replica non-diagonal saddle point.
Let us now consider the calculation of the function D 2 explicitely. First we note, that the saddle point equations eqs. (2.6) are invariant with respect to the replica dependent time shift τ → τ + τ a . This transformation does not change the time dependence in the replica-diagonal elements of saddle point solutions in eq. (2.9), it influences the replica off-diagonal elements though. As we show below, the time-shift transformation is crucial to obtain the correct time dependence for the correlation functions calculated at the replica non-diagonal saddle point. Consider the correlation function between Majorana fermions on different sites (i = j) The correlation function eq. (C.1) in the replica formalism reads where . . . Φ denotes the average over the reparametrization fluctuations, and G ab (τ − τ ) denotes the saddle point replica-off-diagonal Green's function. The time-dependence in eq. (C.2) contradicts the quantum mechanical result, which predicts the dependence of correlation function on the differences τ 1 − τ 2 and τ 3 − τ 4 . However, one can restore the correct quantum mechanical dependence of the correlation function eq. (C.2) by making appropriate time-shifts in each replica. Namely, for each pair of times belonging to the same replica, the time-shift has to be chosen in such a way, that the shifted times are symmetric with respect to zero. Specifically, in eq. (C.2), the times τ 1 , τ 2 , belonging to the replica a, are shifted by c a = −(τ 1 + τ 2 )/2, so that after the shift τ 1 = τ 1 − c a = (τ 1 − τ 2 )/2, and τ 2 = τ 2 − c a = −(τ 1 − τ 2 )/2 = −τ 1 . Correspondingly, the times τ 3 and τ 4 are shifted by c b = −(τ 3 + τ 4 )/2. Let us calculate the correlation function eq. (C.2) assuming the following ordering of times: τ 2 < τ 3 < τ 4 < τ 1 . Using the Liouville quantum mechanical representation for the averaging over the reparametrization fluctuations [46], we obtain One can see, that after the time shifts c a and c b , the correlation function depends on the time differences τ 1 − τ 2 and τ 4 − τ 3 only. For τ 1 = τ 4 = τ , τ 2 = τ 3 = 0 eq. (C.3) reduces to the correlation function calculated in ref. [46], namely Let us now consider the higher powers of the site-nonlocal correlation functions, which we define as a product of 2p quantum mechanical averages To facilitate the transition to the limit case, considered in section 6 we adopt the following ordering of times in eq. (C.4)
To ensure the dependence of eq. (C.7) on the differences of times belonging to the same replica only, we choose the time shifts c k as follows With this choice, the combinations of times entering the exponents in eq. (C.7) become Therefore, the choice of time shifts in eq. (C.9) makes the argument of each exponent in eq. (C.7) to depend only on differences of times in the same replica. Note furthermore, that the time-shifts introduced by eq. (C.8) are replica-local, hence they ensure the quantum mechanically correct time dependence for any choice of contractions.

JHEP09(2019)057 D Details of numerical simulations
In this paper, we use exact diagonalization to investigate the Green's function of SYK model. Majorana fermion operators are represented by γ -matrices, which can be constructed by Clifford algebra [42]. Let γ i to be the representation of the operator χ i . Then one can define the γ (N ) i iteratively. When N = 2, one have γ -matrices as following where σ 1 , σ 2 , σ 3 are the Pauli matrices. Assume we have got γ Diagonalizing the SYK Hamiltonian, one can get all the eigenvalues and eigenstates to construct quantities we want to learn.

E Reparametrization fluctuations around the replica non-diagonal saddle point
In this appendix we show, that the replica non-diagonal saddle point generates coupling between the reparametrizations in different replicas, leaving only a single soft mode, in which all replicas have the same reparametrization. Consider the replica non-diagonal saddle point given by eq. (2.9), and consider the replica off-diagonal part of the softmode action (E.1) which is formulated in terms of the reparametrized Green functions with f a (τ ) being the reparametrization transformation in the replica a. Note, that such form of the Green's function is valid only for times |τ − τ | > 1/J. For shorter times, |τ − τ | 1/J, the Green's function in the model with finite strength of interaction J should approach a free Majorana correlator, G ab free (τ ) = −δ ab sgn(τ )/|τ |. Therefore we restricted the domain of integration in eq. (E.1).
If one further changes the time-integration variables to t i = f (τ i ), defines the field ζ a t = [(f −1 a ) (t)] −1/2 and integrates by parts then the action S 2 [f ] can be cast in the following form (for details see ref. [46]) JHEP09(2019)057 where we introduced for any two functions f 1 and f 2 . Taking into account the symmetries of the Green's function, the polarization operator can be represented in the equivalent form, This expression needs to be found only for times |τ − τ | > 1/J and therefore one can omit the action of time derivative on the sign-function in eq. (2.9). Indeed, the resulting δ-function will bring times t 1,2 infinitely close to each other, but these times are excluded from the integration domain in eq. (E.1). Bearing this remark one obtains with g aa =g and g ab = g for a = b according to eq. (A.1). Then the soft-mode action assumes the form Diagonal matrix elements in this sum produce the Schwarzian action for each function f a (τ ), see e.g. derivation in ref. [46], and below we analyze the terms with a = b. These off-diagonal contributions mainly come from the singularity along the curve C in the plane (τ 1 , τ 2 ) defined by the equation f a (τ 1 ) = f b (τ 2 ). In the close vicinity of C one can introduce the new set of coordinates (τ, s), where τ runs along C and s is the direction which is perpendicular to C, as shown on figure 6. Let e τ ≡ (∂ τ τ 1 (τ, s), ∂ τ τ 2 (τ, s)) and e n ≡ (∂ s τ 1 (τ, s), ∂ s τ 2 (τ, s)) be the corresponding tangential and transverse vector fields to the JHEP09(2019)057 curve C, which in the new coordinates is just a straight line s = 0. According to their definitions we have which means that e τ should be orthogonal to the vector f ab = (f a (τ 1 ), −f b (τ 2 )) T . On other hand e n must be parallel to the same vector f ab . Therefore, the normalized tangential and normal vectors to the line C can be defined by the following relations From here it follows that the parametrization of the coordinate line s in the perpendicular direction e n can be written in the following form where we took into account that f a (τ 1 ) = f b (τ 2 ). Now using the exponential parametrization and substituting eq. (E.11) into the integration kernel in eq. (E.7), we obtain . (E.13) Here we took into consideration that the Jacobian of transformation from (τ 1 , τ 2 ) to (s, τ ) variables is unity, because the basis (e τ , e n ) is orthonormal. The choice of the short-time cutoff ∼ 1/J affects the value of the coupling constant in the action eq. (E.13), it does not change the functional form of the coupling term though. Considering infinitesimally close transformations in different replicas, we can expand the denominator in small difference ϕ(τ ) = φ a (τ 1 (τ )) − φ b (τ 2 (τ )), and after performing the integration over s, we obtain (E.14) For the 2nd term in the expression, ∝ ϕ 2 ab , the contour C can be substituted by the straight line (τ 1 = τ 2 = τ ) if one is interested in the Gaussian order only. As to the 1st part, the line integral

JHEP09(2019)057
by itself has a ϕ 2 -contribution which takes into account the deviation of C from a straight line. As shown below the length reads 16) and therefore the final result for the action in the Gaussian order in fluctuations assumes the form which is the last expression in eq. (3.6).
Let us now derive eq. (E.16). We assume that both phases are small (φ a 1 and φ b 1) and then parametrize C by a variable τ as with fluctuations x, y being small in φ's. Their role is to take into account a deviation of the curve from the straight diagonal. From the relation f a (τ + x) = f b (τ + y) we then have (up to 2nd order in φ), where we took into account that next order terms, e.g f a x 2 φ a e φα x 2 , are of cubic order in φ's. There are many ways to parametrize the same curve and thus equation above does not fix x and y unambiguously. For example one can choose With this choice one further needs to evaluate In the 2nd order in φ's one finds which gives us for the line element When performing the integral over dτ we integrate by parts. Taking into account that f a = e φa we finally arrive at as it was claimed above. It is natural that the correction to the geometric length of the curve C is positive, since it deviates from the straight line.

JHEP09(2019)057
F Transition to the two-level regime In this appendix we give a qualitative explanation as to why the transition from the reparametrization-dominated to the two-level regime occurs at a shorter time-scale for the site off-diagonal correlation functions in comparison to the site-diagonal ones, as observed numerically. Consider first the correlation functions D 2 (τ ) and G 2 ii (τ ) dis . In the two-level regime, with the levels denoted as |0 and |1 , we have D 2 (τ ) = 0|χ i |1 1|χ j |0 0|χ j |1 1|χ i |0 e −2E 1 τ dis = | 0|χ i |1 | 2 | 0|χ j |1 | 2 e −2E 1 τ dis (F.1) Here the energy of the ground state is set to zero hence E 1 denotes the energy gap between the ground state and the first excited state. The numerical data in the two-level regime can be fitted with very high accuracy under the assumptions that energies E 1 and the matrix elements M i = 0|χ i |1 are statistically independent Gaussian distributed quantities. Furthermore, to explain the different time scales for the crossover between the reparametrization dominated and two-level regimes, assume the matrix elements for the operators at different sites to be statistically independent of each other, M i M j dis = 0. Then, for the correlation function D 2 (τ ), we obtain where |M | 2 dis = |M i | 2 dis independently of i. For G 2 ii (τ ) dis we obtain under the same assumptions Under the assumption of Gaussian distributed matrix elements the results eqs. (F.2) and (F.3) differ just by an N -independent factor. The crossover from the reparametization dominated to the two-level regime occurs at the time scale, where the contributions of higher energy levels become suppressed by the corresponding energy exponents ∼ e −Enτ . As a toy model, consider the contribution of the next exited level, which we denote as |2 with the energy E 2 . For the correlation function D 2 (τ ) we now obtain Eq. (F.4) is to be contrasted to the 3-level expression for G 2 ii (τ ) dis G 2 ii (τ ) dis = | 0|χ i |1 | 4 e −2E 1 τ + | 0|χ i |2 | 4 e −2E 2 τ dis + 2 0|χ i |1 1|χ i |0 0|χ i |2 2|χ i |0 e −(E 1 +E 2 )τ dis .
(F.5) Now let us make the following assumptions: matrix elements between the ground state and the state |1 , and between the ground state and the state |2 are statistically independent. Furthermore, the sign of the product 0|χ j |n n|χ i |0 , (i = j) is random, hence the average of such a product over disorder distribution is close to zero. It follows from this assumption, that the terms 0|χ i |1 1|χ j |0 0|χ j |2 2|χ i |0 vanish after the average over disorder, in contrast to 0|χ i |1 1|χ i |0 0|χ i |2 2|χ i |0 , which result in the explicitly positive contribution. Eqs. (F.4), (F.5) become (F.7) Here we denoted M 1 = 0|χ i |1 , M 2 = 0|χ i |2 . Taking into account E 1 < E 2 , on can see that the subleading term e −(E 1 +E 2 )τ dis in eq. (F.7) leads to a slower time decay of the site diagonal correlation function G 2 ii (τ ) dis at relatively short times than the decay of the site off-diagonal correlation function, where the above mentioned term is absent.
Another, complementary point of view on the spread of the correlation functions in the crossover region between the reparametrization-dominated and two-level regime can be gained by considering the effective many particle density of states. Namely, assuming the discrete spectrum of energies E n we conclude that the possible energy factors determining the time decay of the site off-diagonal correlation functions can be only the multiples of the energies E n (in our previous example n = 1, 2), such as 2E 1 , 2E 2 , . . .. This is due to the vanishing disorder averages of the matrix elements 0|χ i |n n|χ j |0 dis for i = j. In contrast, for the site diagonal correlation function, the energy factors are built out of all possible sums of pairs of energies, such as E n + E m , for any two states |n and |m . Considering each energy factor as an effective multi-particle energy level, we conclude, that the many-particle energy spectrum contributing to the site-diagonal correlation function is more dense, that is it has lower many-particle level spacing. Now, the deviations from the reparametrization dominated regime should happen at the times which are of the order of inverse many-particle level spacing. According to the considerations above, those times are larger for the site-diagonal correlation function that for the site off-diagonal one. This would explain qualitatively why the spreading of the site off-diagonal curves [D 2p /(2p − 1)!!] 1/p in figure 4 occurs at earlier times than that for the site-diagonal curves G p 1/p dis in figure 1.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.