Replica-nondiagonal solutions in the SYK model

We study the SYK model in the large $N$ limit beyond the replica-diagonal approximation. First we show that there are exact replica-nondiagonal solutions of the saddle point equations for $q=2$ for any finite replica number $M$. In the interacting $q=4$ case we are able to construct the numerical solutions, which are in one-to-one correspondence to the analytic solutions of the quadratic model. These solutions are singular in the $M \to 0$ limit in both quadratic and quartic interaction cases. The calculations of the on-shell action at finite integer $M$ show that the nondiagonal replica-symmetric saddles are subleading in both quadratic and quartic cases. We also study replica-nondiagonal solutions of the SYK in the strong coupling limit. For arbitrary $q$ we show that besides the usual solutions of the replica-diagonal saddle point equations in the conformal limit, there are also replica-nondiagonal solutions for any value of $M$ (including zero). The specific configurations that we study, have factorized time and replica dependencies. The corresponding saddle point equations are separable at strong coupling, and can be solved using the Parisi ansatz from spin glass theory. We construct the solutions which correspond to the replica-symmetric case and to one-step replica symmetry breaking. We compute the regularized free energy on these solutions in the limit of zero replicas. It is observed that there are nondiagonal solutions with the regularized free energy lower than that of the standard diagonal conformal solution.

The Sachdev-Ye-Kitaev model [1][2][3][4] is a quantum mechanical model of N Majorana fermions with disordered interactions that is solvable at large N in the strong coupling limit. It was proposed [2] as a solvable toy model for holographic description of quantum gravity in the AdS 2 spacetime (see [5] for a review). This idea is justified by the fact that the SYK model displays emergent approximate conformal symmetry in the strong coupling regime [2,3,6] and that it exhibits maximal quantum chaos [2,3,7] at strong coupling. The Goldstone mode corresponding to the conformal symmetry is connected to the gravitational mode in the effective description of the Jackiw-Teitelboim gravity in the AdS 2 bulk [3,4,[8][9][10][11]. Its dynamics in the leading order in inverse coupling turned out to be completely solvable [10,[12][13][14], and in the leading order in 1 N all correlation functions of operators dual to the matter fields in the bulk were computed as well [15]. While the precise bulk dual theory is still unknown, the SYK model has already allowed to obtain significant insight in the physics of black holes and wormholes [16,17]. The defining characteristic of the SYK model is the quenched disorder which randomizes the couplings between sites. Under the assumption that the system is selfaveraging, one can perform averaging over the disorder by introducing replicas, and obtain the path integral in terms of auxiliary fields G αβ (τ 1 , τ 2 ) and Σ αβ (τ 1 , τ 2 ), where α, β are the replica indices [2,4,12]. The results for SYK regarding the large N solution and thermodynamics, are obtained under the assumption that the auxiliary fields are diagonal in replicas. This assumption is justified by the exact diagonalization numerics [3,12,16,18], and by a physical qualitative argument [7,19,20] that prohibits realization of replica-nondiagonal behavior. Besides that, the work which studied spin glass phases (which are usually realized as particular replica-nondiagonal saddle points of the free energy path integral) [1,[19][20][21][22][23] found no numerical evidence of glassy behavior in the fermionic SYK model and provided several analytic arguments of why there shouldn'be any. However, all these considerations do not exclude the existence of the replica-nondiagonal saddle points of the path integral at either finite or zero replicas (in the case of the free energy). Meanwhile, recent work [16,24,25] hints that replicanondiagonal saddle points in annealed quantities, such as the spectral form factor, are responsible for manifestations of quantum chaotic behavior in black holes and possible recovery of information from the black hole, which is lost at the semiclassical level, via holographic duality.
Motivated by these points, in the present work we study the replica-nondiagonal large N saddle points of of the SYK model at general replica number. We start off with the q = 2 variant of the model (where q is the degree of the interacting Hamiltonian). We obtain a family of analytic replica-nondiagonal solutions of the saddle point equations of the disorder-averaged partition function for M replicas of the SYK chain, where M can be understood as an arbitrary real non-negative number. These solutions have an important property of being singular at M → 0. This means that these saddles do not contribute to the free energy. Computing the on-shell action on these solutions at finite replica number, we show that for analytically continued 0 < M < 1 there is a nontrivial phase structure of the path integral, however for M > 1 the standard diagonal solution always dominates. As a next step, we use these analytic solutions in the q = 2 model to construct numerical solutions to the exact saddle point equations in the interacting q = 4 variant of the model. These numerical solutions also exhibit the singular behavior at M → 0 and J = 0. In the q = 4 case the replica-nondiagonal saddles are also subleading in the replica partition function at M > 1.
In the second part of the work, we focus on the strong coupling, or IR limit of the SYK model. In this limit one can find analytic replica-nondaigonal solutions in the interacting model either at finite replica number or in the limit of zero replicas. We study the class of solutions, for which the time dependence and replica dependence are factorized. In this case the saddle point equations separate in the strong coupling limit. The solution is then constructed by solving the equation for the temporal part in the same way as in the replica-diagonal case, and by solving an algebraic equation for the replica part. We perform the latter by using the Parisi ansatz [26,27]. The algebraic equation for the replica dependence in the replica limit transforms into an integral equation. To solve it, we restrict ourselves to the step-function ansatz, which corresponds to the one-step replica symmetry breaking. We study the solutions at M = 0 and compute the regularized free energy on the corresponding saddle points.
The paper is organized as follows. In the section 2 we briefly review the main features of the replica-diagonal saddle-point of SYK, and discuss the non-perturbativity of the exact replica-nondiagonal saddles. We start the study of the replica-nondiagonal solutions in the section 3, which is devoted to the quadratic variation of the model. The exact nondiagonal solutions and their properties are discussed. The next section 4 contains a description and the results of the numerical study of exact saddles in the interacting q = 4 model, and also a general remark on the large replica number limit. In the section 5 we switch gears to the study of replica-nondiagonal solutions in the strong coupling limit. Assuming the factorized ansatz, we derive the reduced saddle point equations and explain in detail the general strategy for constructing the solutions and computing the regularized free energy in the zero replicas limit. Subsequently, in the section 6 we construct the solutions in the one-step replica symmetry breaking ansatz and compute the leading contribution to the regularized free energy in the strong coupling limit. In the next section 7 we make some comments about generating other solutions using the reparametrization symmetry and possible holographic interpretation in particular cases. We discuss our results and unanswered questions in the section 8. The appendix A provides a brief introduction to Parisi matrices and derivations of a few formulas used in the main text, and the appendix B contains a few general formulae regarding the computation of on-shell action in different cases.

Setup
The object of our study is the Sachdev-Ye-Kitaev model [2][3][4], which is a theory of N ≫ 1 interacting Majorana fermions in 0 + 1 dimensions 1 . The Hamiltonian is given by Here ψ i are the Majorana fermions, and j i 1 ...iq are totally antisymmetric couplings randomized via the Gaussian distribution: (2.2) 1 In this paper we work in the Euclidean time, unless mentioned otherwise.
To calculate a physical quantity in this model, one has to average over the disorder using the rules which follow from the distribution (2.2): To study physically meaningful quantities, one usually has to average over all realizations of the disorder. The free energy of the model with quenched disorder is given by where Z = Tr e −βH . To simplify evaluating the disorder average, one employs the replica trick, which we write in the form: In the case of integer M on the right hand side there is a path integral over M copies of SYK, so the fields now carry an additional index α = 1, . . . , M . We study the case of integer values of M separately. The limit M → 0 requires an analytic continuation, which will be considered on a particular ansatz. One can calculate the disorder average for the replica partition function and rewrite it in terms of the path integral over O(N )-invariant auxiliary fields with replica indices. The derivation was presented in detail in [4,12]. After this procedure, one obtains the following expression for the replica partition function: where G is the bilocal field which has the meaning of the Majorana fermion two-point function, and Σ is the auxiliary bilocal field 2 which has the meaning of the fermion self-energy. These bilocal fields are supposed to satisfy the antisymmetry condition In the present work we essentially study the saddle points of (2.6) for different values of M . The saddle points of the path integral are defined by the following equations:

Review of the replica-diagonal solution
The replica partition function (2.6) has a family of replica-diagonal saddle points, which have been extensively studied in the literature [2-4, 7, 12, 15]. One assumes the ansatz The quenched average in this case coincides with the annealed average up to subleading orders in 1 N expansion [4]: 11) which means that up to subleading orders in 1 N one can take off the replica indices in the path integral (2.6): (2.12) At large N the asymptotic of the RHS of (2.12) is given by saddle points contributions where the saddle points of (2.12) are given by the Schwinger-Dyson equations for melonic diagrams [2,3] 1 ∂ τ −Σ =Ĝ ; (2.14) where in the first equation hats denote the integral operators with the kernels defined by the corresponding bilocal fields. The solutions of these equations in general case can be constructed numerically, which was done in the previous work [3,4,12,16], and analytic solution is known in the IR/strong coupling limit βJ ≫ 1.

Strong coupling limit
Substituting (2.15) into (2.14) and taking ∂ τ → 0, one obtains the equation (2.16) The solution for G of the equation (2.16) has the form of the conformal propagator on the circle: where ∆ = 1 q is the conformal dimension of the Majorana fermion and In the frequency space the conformal propagator on a circle has a form [8,28]: where ω n are Matsubara frequencies: At zero temperature (2.17) and (2.19) reduce to where b is the dimensional constant fixed from (2.18). The equations (2.14), (2.15) are invariant under time reparametrizations τ → f (τ ) in the strong coupling limit, provided the bilocal fields transform as follows (we assume without loss of generality that f is a monotonically increasing function): Acting with these transformations on the solution (2.17), one can generate the full infinite-dimensional manifold diff(S 1 ) SL(2, R) of the replica-diagonal saddle points.

On-shell action for replica-diagonal solution
The replica-diagonal on-shell action is given by where G and Σ solve (2.14) and (2.15). The conformal limit is obtained by neglecting the time derivative. The s 1 and s 2 can be simplified and we can rewrite them in the equivalent forms In the conformal limit both these pieces in the on-shell action diverge and have to be regularized. It is expected that the renormalization can be performed pertrurbatively (in 1 βJ) by reinstating the time derivative in s 1 and evaluating the corresponding counterterms, and the resulting renormalized action equals to the on-shell action, evaluated on the solution of exact saddle point equations [3].

Non-perturbative nature of the replica-nondiagonal correlators
Here we will show that one cannot obtain a replica-nondiagonal large N solution in perturbation theory over the free fermionic theory. We will work in the frequency space, and we also fix q = 4. In the free case we have J = 0, and the saddle point equations (2.8)-(2.9) are solved by We introduce the dimensionless parameter λ and look for the solution by perturbing the free UV fixed point: αβ (ω) + . . . , . (2.32) The field Σ in this case is expanded as follows: αβ )(ω)δ αβ (2.33) αβ )(ω)δ αβ + J 2 λ(g (1) αβ * g (1) where the replica matrices are always multiplied a-la Hadamard, i.e. component-wise.
Note that the leading possible replica-nondiagonal contribution to Σ is of order λ q−1 . Substituting all this into the saddle point equation (2.8) and equating the powers of λ, we arrive at an infinite system of linear inhomogeneous integral equations for g (k) . For example, for λ 1 we obtain It is clear that the solution for g (1) of this equation is replica-diagonal. Having solved this equation, one can substitute the solution into the λ 2 equation, which would then allow to solve for g (2) . However, since everything in the equation will be replicadiagonal, the solution for g (2) will also be replica-diagonal. Using this expansion one can construct the solution up to any finite order in λ, and it will thus always remain replica-diagonal, ultimately because the free fixed point G f is replica-diagonal. This means that we cannot obtain a replica-nondiagonal solution as a low-energy limit of a replica-diagonal solution, and any replica-nondiagonal large N solution would be non-perturbative in (βJ) −1 .

Nondiagonal saddles in the q = 2 model
In this section we consider the q = 2 variant of the SYK model. We present a class of simple exact replica-nondiagonal solutions and study their properties.

The solutions
Let us consider the q = 2 case and show that there are replica-nondiagonal solutions of the saddle point equation. The saddle point equations (2.8)-(2.9) in the q = 2 case condense to a single SD equation, which in terms of replica matrices is written as where ω is the Matsubara frequency, and I is the unit matrix in the replica space. We assume the replica-symmetric ansatz The equation (3.1) turns into the pair of equations for G 0 and G 1 : In order to obtain the fermionic solutions, we have to impose the antisymmetry condition G αβ (ω) = −G βα (−ω) , (3.5) which in terms of the replica-symmetric ansatz simply means that G 0 and G 1 have to be odd functions in frequency and time domains. The equations (3.3)-(3.4) are readily solved. There are two replica-diagonal solutions: and two replica-nondiagonal solutions Let us make a few remarks about these solutions. The first solution coincides with the solution presented in [3]. It admits the J = 0 limit, where its leading asymptotic is given by the free correlator 1 −iω . The second solution is the one, which gives subleading saddles, discussed by Cotler et al in [16]. All solutions are pure imaginary and There are also the relations In regards to the replica-nondiagonal third and fourth solutions, the important property worth pointing out here is that they are singular in the M → 0 limit. Also, one can check that these solutions are singular in the free limit J → 0, which confirms the nonperturbative nature of these solutions, discussed in the sec. 2.2. We will see that these properties remain for the class of numerical q = 4 solutions that we studied. We conclude this subsection with a comment about other solutions of the equation (3.1). This is a quadratic matrix equation, which means that in principle one can find and classify all of the solutions at any fixed M . Their general form can be found by rewriting the equation (3.1) as follows 3 : The general solution will have a form where X is a matrix which parametrizes a particular solution, such that: There is no explicit dependence on M in (3.18), however the diagonal component of equation (3.19) has a sum of M terms equal to 1. That means that the individual nontrivial components of X should contain the 1 M dependence to compensate. Because of this argument, we expect that every solution will have singularity at M → 0.

On-shell action
Now we turn to the study of contributions of the nondiagonal saddle points, described above, to the replica partition function. Since the saddle point equations for different frequency modes decouple, we can consider the density of the action, which we denote by ρ and define at zero temperature as follows: Here V is the regularized volume. At finite temperature the definition is generalized by setting V = 2π and ∫ dω → ∑ ωn . The decoupling of the saddle point equations means that for every allowed frequency one can in principle choose any of the four solutions obtained above. The question, in which we are interested here in this section,     is whether there are any saddle points that would dominate over the replica diagonal saddle. Because of the frequency decoupling, to check this fact it is enough to compare the action density ρ evaluated on different roots. Using the formulae (B.6),(B.8) at q = 2, the action density is written as (3.21) We denote the contributions from different solutions (3.6)-(3.12) as If we see an interval where ρ (3) < ρ (1) (ρ (4) < ρ (1) ), then that means that there are solutions that have action lower than the replica-diagonal. One such solution can be constructed by selecting the 3-rd (4-th) root for the frequencies from the interval where the above inequality is true, and choosing the 1-st root for the rest of the frequencies.
Thus we focus on the study of the function ρ(ω) in the remainder of this subsection. Let us first estimate the contributions of the first two diagonal solutions to the on-shell action. For these roots l 2 = 0 and: We can expand in terms of small ω > 0. In this case we get Let us now present the contributions to the action density the last two solutions.
For small frequencies ω → 0 (ω > 0), we get Comparing (3.30) with (3.27) we see in the leading order the asymptotic coincides, except for phase contribution, which would be proportional to 2πi in the total action and thus inconsequential. However, in the subleading ω 1 order and higher there is difference. As hinted by this asymptotic and confirmed by the plots of the exact expressions on Fig.1, the M = 1 is a threshold value which distinguishes between two different types of behavior of saddle points:  (4) , and for M < 2 we have ρ (3) > ρ (4) . This is also illustrated on the plots Fig.2D-F we plot the difference between the contributions of two nondiagonal solutions ρ (3) − ρ (4) as a function of ω and M . On Fig.2A-C we plot the difference ρ (1) − ρ(non-diag) 4 . Besides the observations mentioned above, from these density plots it is evident that the IR region seems to be more robust in the singular M → 0 limit, rather than the UV. One can interpret this as a hint towards the fact that the singular behavior in M → 0 limit of nondiagonal solutions is the UV effect in the SYK model. We explain more evidence for this in other sections of the paper.

Exact nondiagonal saddles in q = SYK: numerical study
Having found nondiagonal solutions in the q = 2 model, we now turn to study the interacting q = 4 model. In this case the saddle point equations (2.8)-(2.9) cannot be solved analytically in general, so we construct the solutions numerically. As was done in the previous section, we assume the replica-symmetric ansatz (3.2). In terms of independent variables the saddle point equations read We solve the equations numerically at finite temperature, with β = 2π.

Comments on the method
We solve the system of integral equations (4.1)-(4.3) by iterating them. We use the approach employed in [16] in studies of subleading replica-diagonal saddles. The main idea is to start iterations with q = 2, using a particular solution of the q = 2 model as a trial functions, and gradually increase q from 2 to 4 during the procedure. Let us know discuss the procedure in more detail.
Initial condition. The trial functions for G 0 and G 1 are constructed by choosing one of the four solutions (3.6)-(3.12) for every Matsubara frequency. We want the resulting solution in the interacting model to have the asymptotic behavior in the UV region that would correspond to the free theory, so we only consider the q = 2 trial functions for which ∃n such that G 0,1 (ω n ) = −G 0,1 (−ω n ) = G (1) 0,1 (ω n ) ∀n ≥n. In this case for any n <n we can choose G 0,1 (ω n ) = −G 0,1 (−ω n ) to be equal to any of the four solutions (3.6)-(3.12).
Iteration procedure. We divide the iterations into two stages.
1. We start first stage of iterations at q = 2. At each iteration, q is increased by some small amount. At every step Σ 0 and Σ 1 are computed in the position space using fast Fourier transform for G and equation (4.3), and then the inverse fast Fourier transform is performed on Σ 0,1 . Then G 0 (ω) and G 1 (ω) are updated according to the weighted rule (as also used in [3]) of the form where 0 < x < 1 is the weighting coefficient andG is defined by solving the equations (4.1)-(4.2) in terms of G 0 and G 1 : ; (4.5) .
At this stage we keep the weight fixed, and the procedure is finished when q reaches 4.
2. The second stage of iterations is performed at fixed q = 4. Its purpose is to tune the solutions, obtained in the previous stage, to the desired precision. In this procedure we take the approach of [3] and control the L 2 -norm of the solutions between successive steps decreasing the weight x every time the ∆G 0,1 2 starts increasing. The procedure is completed once ∆G 0,1 2 reaches zero (up to desired numerical accuracy).
The main limitation of our approach is that only real-valued stable numerical solutions can be obtained. This puts limitations on making connections with the analytic solutions in the conformal limit, which we discuss in the section 5) and thereafter.

The results
When studying the saddle points of the replica partition function at finite M , the obtained solutions indicate that for every q = 2 solution, that we choose as initial condition as discussed above, there exists a solution of the q = 4 model. This generalizes the observation, made in [16] for the replica-diagonal solutions, to the replica-nondiagonal symmetric case. The solutions shown on Fig.3 are obtained by iterating from the following q = 2 solutions: for all other n.
0,1 (ω n ) for all other n.  We have also studied this class of solutions in the limit M → 0. For this purpose we add a third stage of iterations, where we keep q fixed, but change the value of M from some finite initial value to zero during iterations. Evaluating ∆G 0,1 2 and the left hand side of the equations of motion, we observe that the sequence of functions obtained this way fails to converge to any solution of the saddle point equations, other than the standard replica-diagonal saddle. When starting from a replica-nondiagonal solution at finite M , we observe that the iterated functions develop discontinuities as M → 0. Thus, from our numerical evidence we conclude that the singular behavior in the M → 0 limit, that we see in the analytic q = 2 solutions, persists for nondiagonal solutions in the q = 4 model. Therefore, the main result of our numerical investigation is that at finite replica number we get an infinite number of nontrivial replica-nondiagonal saddle points, whereas in zero replicas limit we do not obtain any replica-nondiagonal solutions of the exact saddle point equations.

On-shell action
We compute the on-shell action on the replica-nondiagonal solutions at finite M defined by the formula (B.1). it appears that in the q = 4 case all the nondiagonal saddles that we have constructed for M > 1 are subleading, similarly to the q = 2 case. We study the difference where S(standard) is the value of the on-shell action on the standard replica-diagonal saddle (times M ), and S(nondiagonal) is the value of the action on a particular replicanondiagonal solutions. On Fig.4 we plot 2 N ∆S as a function of the label n of a Matsubara mode, which was chosen to be other than G (1) in the q = 2 trial function of the corresponding numerical solution. For the nondiagonal G (3) and G (4) branches ∆S decays with n according to what appears to be a power law n α with the exponent α determined by M (as well as coupling). We also confirm the approximately linear decay on the diagonal G (2) branch, which was stated in [16].

Remark on the large replica number limit
We conclude our discussion of exact replica-symmetric nondiagonal saddle points in q = 2 and q = 4 models by considering the large replica number limit, M → ∞. Namely, we note that the solutions in this limit become replica-diagonal. In the q = 2 case this is evident from the analytic solutions (3.9)-(3.12). From these formulae we see that the nondiagonal terms G So in principle one can treat the nondiagonal terms as the 1 M -corrections. The diagonality of the solutions in the M → 0 limit leads to the following identity: This can be interpreted as a sort of self-averaging relation for the partition function for large M and N .

Replica-nondiagonal solutions at strong coupling
We have found in the above section that the limit of zero replicas is singular in the replica-nondiagonal solutions of the exact saddle point equations, and thus there are no replica-nondiagonal solutions at M = 0 in the class which we studied. Now we will consider the saddle point equations in the strong coupling (IR or conformal) limit, where we omit the time derivative. By neglecting the UV source term, which is diagonal in replicas, we thus allow for a much bigger set of possible solutions. The aim of the present and the subsequent sections is to construct solutions of the SYK model in the strong coupling limit in M → 0 limit.
In the strong coupling limit βJ ≫ 1 at finite replica number M the saddle point equations (2.8), (2.9) which follow from (2.6) read: The saddle point equations (5.1),(5.2) in the strong coupling limit are invariant under the transformations, which are induced by separate time reparametrization τ → f α (τ ) in every replica [29]: In other words, in the general replica-nondiagonal case the emergent conformal symmetry extends to the group diff(S 1 ) ×M .

Separation of variables in the IR limit
We are going to study the solutions of saddle point equations at strong coupling using the particular ansatz, where the time and replica dependencies are factorized [29]: The main advantage of using the ansatz (5.5) is that we can construct analytic solutions. First, we substitute Σ αβ from the second equation (2.9) into (2.8), and we are left with the equation for G: We substitute the factorized ansatz (5.5): The general solution to (5.7) is given by a matrix P αβ and a function g(τ ) that satisfy the matrix equation and the integral equation here C is an arbitrary non-zero constant. We take the antisymmetric conformal invariant solution of (5.9) where b is defined by (2.18). Due to the requirement G αβ (τ, τ ′ ) = −G βα (τ ′ , τ ) the matrix P is symmetric P αβ = P βα . We note that the equation (5.7) has the scaling symmetry under transformations Thus we have a scaling degree of freedom which can be fixed arbitrarily. There are two convenient ways two impose the scaling symmetry fixing condition.
• Normalization constraint. One can fix C = ∑ β P q βγ = 1. The off-diagonal part of the equation (5.8) reduces to This equation is treated on equal footing with the normalization constraint, and they together determine the matrix P .
• Diagonal constraint. In the present work we instead impose the condition when we fix P αα = 1. In this case one first solves the equation for non-diagonal components of P , which has the form Eq.(5.12), and then computes C to completely determine the solution according to This approach is more convenient for study of specific solutions, and we employ it throughout the paper.
The previous considerations in fact mean that ultimately we deal with the normalized replica matrix, such that The degree (q−1) is then understood as the matrix degree with respect to the Hadamard multiplication, and we write the equation (5.8) as follows: where P ○r = P ○ P ○ ...P r , and I is the identity matrix.
As a side remark, we note that the equation (5.17) with C = 1 can be interpreted as the saddle point equation for the 0-dimensional version of the SYK model [30]. Thus in principle one can treat the factorized ansatz (5.5) as a sort of dimensional reduction of the SYK model.
We also note that the factorized ansatz (5.5) breaks the reparametrization symmetry (5.3)-(5.4) symmetry down to a single copy of diff(S 1 ), which acts in the same way as in the replica-diagonal case, see eqs.
In particular, these transformations can alter time dependence of replica-nondiagonal components and thus lead to physically more interesting solutions. An example is considered in the section 7.1.

The approach
To solve the equations (5.8), we use the Parisi ansatz [26,27] for the matrix P . For the definition and properties of Parisi matrices, see appendix A. The Parisi matrix P of the rank l is characterized by the parameters a 0 , . . . , a l . The key properties of the Parisi matrix form that make it especially suitable for solving the equation (5.8), are the following: 1. A Parisi matrix satisfies the constraint β P q αβ = C ∀α (5.21) identically, for any q with some C. This is easy to see because every line and every column of a Parisi matrix contains all of its parameters, so the sum across every line and column is the same.
2. The Parisi matrices form an algebra with respect to both direct matrix product and the Hadamard matrix product (the proof is presented in the section A.3). This guarantees the consistency of the matrix equation (5.17) (see representation (B.12)).
For different possible configurations, we use the terminology analogous to the context of spin glass solutions: Our aim is to obtain and study solutions in the limit M → 0. We approach the problem of finding the solutions 5 . We will take the limit M → 0 first, directly in the saddle point equation (5.17) and in the on-shell action. Then we will find the solutions and calculate the (regularized) on-shell action on them.
In general, the strong coupling limit of SYK contains UV divergences, and they will appear in the action, evaluated on an IR solution. Therefore, if one wants to understand the role of solutions in this limit when submerged into complete SYK model, one has to perform the renormalization, by fitting to the numerical solutions of the exact equations 6 . Because we were not able to find solutions of the exact saddle point equations in zero replicas limit as discussed in section 4, we expect that for any replica-nondiagonal solution in the IR limit this renormalization will bring the value of the on-shell action to be equal to the standard saddle value. Thus we do not expect that the nondiagonal solutions that we construct in M = 0 case in the strong coupling limit make any contribution to the physical free energy in the complete SYK model. Nevertheless, we will be interested in calculating the regularized free energy in the strong coupling limit on nondiagonal solutions, and study its difference from the regularized free energy on the diagonal conformal solution in order to analyze the dominance of saddles in the leading order of the strong coupling expansion.
We employ the strategy, outlined above, in the one-step replica symmetry breaking case and discuss the corresponding solutions in detail in the next section 6. The remainder of this section is focused on technical aspects of regularization and calculation of the on-shell action and of the limit M → 0.

Pfaffian factorization
To obtain the expression of on-shell action for factorized solutions from the partition function (2.6), we first need to evaluate the Pfaffian term. In the IR limit ∂ τ → 0 the factorized ansatz for replica-nondiagonal solutions allows for the factorization of the Pfaffian. To prove it, we can write We can diagonalize the Σ αβ (τ, τ ′ ) by making the transition to the frequency space. We write where ω n are Matsubara frequencies defined in (2.20). Then the quadratic form reads On the factorized solution (5.5), we write where Σ c (τ, τ ′ ) = J 2 g c (τ, τ ′ ) q−1 , and g c is the conformal propagator given by (2.17). Substituting this and (5.27) into (5.24), we evaluate the path integral as The factor depending onΣ is the same as in the replica-diagonal case. The second factor is the contribution from replicas, it is an infinite degree of the determinant of a Parisi matrix. The latter is calculated in the M → 0 limit in section (A.4), see eq.(A.43). To calculate its contribution in the action, we only need to introduce an appropriate regularization by introducing a cutoff at some large n in the product.

Regularization
To finally separate out the contribution of the replica matrix in the Pfaffian, we need to regularize it. In the present work we use a direct regularization, which is a hard cutoff in the frequency space such that the validity of the strong coupling regime is preserved: and The restriction (5.31) also supports the validity of the strong coupling, or IR, regime. Indeed, the contribution of the free propagator (the ω-term) in the denominator of the LHS of (2.14) is suppressed as compared to Σ(ω) that at large β can be estimated as and ω n < Σ(ω n ) corresponds to the bound (5.31).
Another regularization was discussed in [28], it is essentially equivalent to the exponential cutoff in frequencies.

On-shell action on factorized solutions
We start with the on-shell action for the partition function (2.6) at finite M in the strong coupling limit, Substituting for the solution we get Using the factorization property of the det and the identity α,β P ○q αβ = CM , we obtain where s 1,c and s 2,c are defined in (2.28) and (2.29) respectively. Note that we use the same regularization in the polynomial term (which amounts to (5.35), as discussed above for the Pfaffian term. The new term unique to replica-nondiagonal factorized solutions is introduced: Note that it is proportional to the UV cutoff parameter d f , which is defined in (5.35).

5.5
The M → 0 limit

Replica symmetric case
Let us first consider the simplest case of a Parisi matrix to illustrate our approach for taking the zero replicas limit. We consider the Parisi matrix P of the first level (as explained in appendix A), or the replica symmetric ansatz. Let us take P αα = a, α = 1, 2, ..., M and P αβ = A if α ≠ β. Here a and A are two generically complex-valued numbers. We have The limit M → 0 in this case is taken straightforwardly:  and in the limit M → 0 take the form Fixing the scaling freedom by setting a = 1,we get the final equations for the m = 0 case: For q = 2(n − 1) the equation (5.50) has n pairs of conjugated roots (b (i) ,b (i) ). Each such pair gives the contribution to the regularized free energy, determined by the real part of s 3 , see discussion below in sec. 5.5.4 for more details. Taking into account that the cutoff is d f ≃ βJ, we can compute the normalized difference between the regularized free energies of the replica-diagonal and nondiagonal solutions: (1)     Fig.5A, and the absolute value as a function of q is plotted on Fig.5B. The values of f q for other replica-symmetric solutions are plotted on Fig.6A, and the dependence of the f q for the lowest replicasymmetric saddle on q is presented on Fig.6B.

General Parisi ansatz at M → 0
To take the limit M → 0 for more general structure of the Parisi matrix, we use the standard trick of [26]: we introduce the Parisi function of a continuous variable instead of the Parisi matrix parameters, thus performing the analytic continuation from integer values of M to arbitrary positive real number: a j → a(u) Taking the limit M → 0, we write we arrive at the following equations: where u ∈ [0, 1]. Finally, we fix the scaling freedom by imposing P αα = a 0 = 1, or, equivalently, dividing (5.65) and (5.66) by a q 0 . We arrive at Thus, these integral equations are the final form of the matrix equation (5.17) in the zero replicas limit.

A comment on the q = 2 case
For q = 2 the equations (5.68),(5.67) have the form Taking the derivative on the first equation, we get Therefore the replica-diagonal solution is the only valid solution of the q = 2 SYK model in the zero replicas limit. We note that the proof works only for smooth functions a(u).
For the discontinuous function a(u), we have checked that this also is the case on the step function ansatz, and expect this to be true for any function.

Contribution to the regularized free energy
The free energy is expressed from (2.4) using (2.5): To calculate it in the large-N approximation, we have to find the saddle point configuration of Z M ∼ exp(−S M ) in the space of replica bilocal fields G αβ and Σ αβ with the minimal value of the real part of the on-shell action. On the factorized solutions, we can have one or several saddle points with the same value of real part of the action. We are specifically interested in the free energy of factorized replica-nondiagonal Parisi saddle points in the conformal limit. The regularized on-shell action in this case, as dictated by (5.41), separates into the contribution identical to the replica-diagonal on-shell action in the conformal limit plus a contribution from the Parisi replica matrix: where we use the cutoff consistency assumption that d f = βJ π, and the matrix Q on the replica space is defined as The matrix Q has the Parisi form and is expanded as (A.12): where the Parisi parameters are defined as follows: In the continuum representation the matrix Q is defined by the constant q 0 and the Parisi function q(u), which are defined as (with a 0 = 1 taken into account) To compute (5.76) in the replica limit M → 0, we use the formula for the Parisi matrix tracelog in the continuum representation [26]: Substituting (5.77) into (5.82), we obtain the expression: With C substituted using the equation (5.68), this formula establishes the resulting general expression for the contribution of the replica structure to the on-shell action on a particular solution for the Parisi function a(u) 7 .
To be able to compute the free energy, what is left is to describe which saddle points actually contribute in the path integral. We find them by solving the equations of motion for the Parisi matrix P . The subtlety here is that we can have multiple saddle points with equal absolute value of the integrands. Each saddle point gives a contribution to the free energy which reads where k labels the saddle point. If we have a family of n saddle points with free energy values F saddle (k) such that Re F saddle (i) = Re F saddle (j) ∀i, j = 1, . . . , n , (5.87) we have to sum over them to obtain the full expression for the free energy: Note that in the case of integer M the l.h.s. of the equation (5.17) is a polynomial with real coefficients. Therefore is P 1 is a complex-valued solution of the equation, then the matrix P 2 = P * 1 is also a solution, as we saw in the replica-symmetric particular case in sec. 5.5.1. Therefore for every complex saddle point we will also have the complex conjugated saddle point contributing in the sum in (5.88). Focusing on the case of two complex-conjugated saddle points in the sum, the formula for the free energy is written as Recall that the saddle point value of the free energy is determined by the on-shell action S M in the M → 0 limit, see Eq.(5.86). The S M is expressed by the formula (5.41). First of all, let us note that the on-shell action is proportional to N . The first term in the formula (5.89) is a real part of that action, and is also therefore an extensive contribution to the free energy. However, the second term is not proportional to N .
In the large N limit the second term can be neglected. Therefore, we are left with the real part of the on-shell action (5.41) defining the free energy at large N . The s 3 piece is responsible for non-trivial contributions to the free energy of the replica-nondiagonal factorized solutions. Hence on specific solutions we are most interested in the quantity Re s 3 , (5.90) and, in particular, in the sign of ∆F describing the shift of the replica non-diagonal solutions in respect to the diagonal one. The negativity of ∆F for given pair of solutions means that the replica non-diagonal solution is the dominant one. This is the quantity in the rigid strong coupling limit where the replica matrix P in the factorized ansatz introduces discrepancy with the replica-diagonal result.

One-step RSB solution
Having established the formalism and derived general expressions for equations of motion and on-shell action in the previous section, we are now ready to study solutions more specifically. In this section we focus on the SYK model with q = 4. We study the solutions of the equations (5.68),(5.67). We restrict ourselves to the solutions for a(u), which can be described by the one-step replica symmetry breaking ansatz, in analogy with the one-step solutions in the spin glass systems [26,27], such as the Sachdev-Ye model [19]: In this formula µ is a free parameter, to which we will refer as the breakpoint (again, in analogy with [19]). The moments of a(u) and integrals which contribute to (5.68), (5.67) are evaluated on the ansatz (6.1) as follows: The equation (5.68) reduces to Meanwhile, the equation of motion (5.67) decays into two equations which we get separating the coefficient in front of the step function: We take µ as the free parameter and solve these equations for A 0 and A 1 . Equations (6.7)-(6.8) are algebraic equations of the 4-th order, so we have 16 solutions. We want to pick solutions that describe saddle points of the model, i.e. ones that also solve the equation (5.7) in the M → 0 limit. We compute the contribution to the regularized free energy (5.90) on these solutions For the tracelog we use the formula (5.85), which on the one-step ansatz assumes the following form:  • All of the complex saddle points, which give the same real part of the action, are organized in mutually conjugated pairs, and we didn't find any instances of multiple pairs of saddle points giving the same real part of the on-shell action, so the formula (5.90) is applicable.
• While the general formula (5.89) does not necessarily guarantees that the regularized free energy is real-valued, we found this to be the case for all solutions we studied, except for one which we mention below.
The constant solutions correspond to one step function with A 1 = 0, see (6.1). We have two kinds of solutions in this class: 1. Replica-diagonal (paramagnetic) solution: A 0 = 0, A 1 = 0. 2. Replica-symmetric complex-valued solutions. This case is reduced to considerations performed sec. 5.5.1, see table 1.
Replica symmetry breaking solutions. There are also two kinds of solutions with A 1 ≠ 0: 1. A 0 = 0. In this case A 1 is a solution of the equation which follows from (6.8): This equation has one real and two complex mutually conjugated solutions: ; (6.13)  where The regularized free energy for the real solution labeled by (1) is complex valued. The location of solution (6.13) on the complex plane and its contribution to the free energy is presented on Fig.7. Note that the free energy diverges at the real solution near the limiting case µ = 0.
The solutions (6.14) are mutually complex conjugated. Their trajectories on the complex plane parametrized by µ and their contribution to the regularized free energy is presented on Fig.8. We see that at µ = 0 the positions of these solutions coincide with the positions of pure constant complex solutions (see table 1 at q = 4), shown by green points. This is not surprising since for µ = 0 our step function solution in fact corresponds to the constant solution.
2. Second group of RSB solutions is characterized by non-zero both A 0 and A 1 . In this case the degree of the algebraic equations (6.7),(6.8) is too high to obtain explicit analytic solutions, but they can be solved numerically. We plot these solutions on the complex plane on Fig.9. Among these RSB solutions there are some that turn into the real replica-symmetric solution with A 0 = a 0 at either µ = 0 or µ = 1, whereas others reduce to the complex replica-symmetric solutions. We have calculated the contribution to the regularized free energy of the remaining solutions. The dependence of ∆F on the breakpoint parameter µ is presented on Fig.10. We see that there are three local minima below the replica-diagonal value.  For example, we can make the time dependence of replica-offdiagonal components of the field G αβ differ from the time dependence of the diagonal components if we act with the set of diffeomorphisms V α ∈ diff(S 1 ) which is arranged in such a way that In this case all diagonal components will have the same time dependence, whereas the time dependence of the off-diagonal components will be different, because of the extra SL(2, R) transformation which acts only on one of the times.
Example. Let us fix M = 2. We will work in terms of dimensionless times on the thermal circle Consider the following pair of diffeomorphisms, where the first transformation is identity and the second one is a fraction linear transformation on the circle: In this case the diagonal components G 11 and G 22 are invariant under these transformations and are given by However, the off-diagonal part transforms non-trivially: In order to analyze the dependence on real time, we set θ 1 = 0 and e iθ 2 = e 2π β t 2 . Then, we obtain This is the decay law of the off-diagonal component. The exponent does not dominate until the characteristic time scale of Since the SL(2, R) is non-compact, by fine tuning the parameters of the Möbius transformation, we can make this time scale as long as we like, so that For comparison, the diagonal components behave like (7.10) and it starts decaying right away. This means that at short times t 2 < t * 2 the transformed solution is frozen, behaving like a kind of spin glass. However, after t * 2 the thermal fluctuations destroy the spin-glass-like configuration with the same speed as a regular replica-diagonal configuration. The emergent (extended) conformal invariance suggests that the factorized solutions might have a holographic interpretation in terms of nearly AdS 2 gravity, like the replicadiagonal solution does [9,11,24]. The purpose of this section is to explore this correspondence on a toy example of factorized solution for M = 2.
In this case solving the equation (5.8) for the q = 4 case gives a solution with P 11 = P 22 = 1 ; P 12 = P 21 = i ; (7.11) We now apply a pair of reparametrizations to the solution using (5. 3): Each of these transformations belongs to the unbroken SL(2, R), so diagonal elements of G are unchanged, whereas the off-diagonal elements transform non-trivially. The transformed solution is given by 2∆ ; (7.14) When performing the analytic continuation to the Lorentzian signature, it is evident that (7.14), (7.15) are, correspondingly, one-sided and two-sided correlators of a Majorana fermion in the thermofield double state [24,34,35]. The transformations (7.13) make it explicit that the two-sided (off-diagonal) correlator can be obtained from a onesided (diagonal) correlator by moving one of the endpoints halfway along the thermal circle [35]. Thus, a replica-nondiagonal large N saddle point of the SYK model with M = 2 replicas in the conformal limit can describe the purification of a single replica of SYK model at finite temperature in the form of the thermofield double state.

Comparison with semiclassical holographic computation
We can compare the factorized solution (7.14), (7.15) with what we can get from a holographic computation of correlators in the AdS 2 spacetime. We consider the Lorentzian AdS 2 spacetime, which is described in embedding space of signature (− − +) as a hyperboloid defined by the equation [9,24] − T 2 1 − T 2 2 + X 2 = −1 . (7.16) We are interested in the Schwarzschild coordinate patch in the AdS 2 . Let us denote by symbols L and R two boundaries of the AdS 2 . In this case there are two corresponding choices of the Schwarzshild patch, which are described by the parametrizations [24]: Here r, t ∈ R. The horizon radius is related to the temperature as R = π β . The induced metric in both cases is the same, (7.20) Suppose we want to calculate a correlator of heavy (scalar) operator O of the dimension ∆ O inserted on the same boundary. In the limit ∆ O ≫ 1, the correlator is determined 8 by the length L of the geodesic anchored on the boundary endpoints: Thus we need the boundary-to-boundary geodesics. The one-sided geodesics give complex-valued lengths because the interval is timelike, but we can circumvent this by making analytic continuation to the Euclidean signature τ = it in the parametrization and T 2 → iT 2 in the embedding space. To find the geodesic length between the points 1 and 2, we can use the relation where ⃗ Y denotes a point in the embedding space, and angular brackets denote the scalar product in the embedding space. We choose the endpoints on the boundary, i. e. so that r 1 = r 2 = r 0 → ∞. Taking the asymptotic with respect to r 0 and subtracting the divergent part, we obtain for one-sided correlator and for the two-sided correlator We see that the time dependence of these semiclassical correlators is captured properly by the 2-replica result (7.14), (7.15). However, it is clear that the structure of the Parisi matrix P will not be captured by the bulk theory on the leading semiclassical level.
In conclusion to this remark we would like to note that in the holographic derivation we performed the analytic continuation to the Euclidean signature while trating it purely formally. It seems plausible that the actual bulk spacetime which would have to correspond to our analytically continued Lorentzian AdS 2 would be similar to the double cone construction, described by authors of [25].

Discussion
In this paper we have discussed replica-nondiagonal saddle points of the replica partition function (2.6) of the SYK model. First, we have studied the exact nondiagonal saddles, analytically in the q = 2 model and numerically in the q = 4 model. Let us summarize and comment on our main results from this study.
1. For q = 2 we found exact analytic replica nondiagonal solutions of the saddle point equations, given by the formulae (3.9)-(3.12). As we have discussed, these solutions are singular in the limit M → 0. It is easy to check that the on-shell action in the zero replica limit is also singular, which means that there is no well-defined limit of zero replicas on this class of solutions. Another significant property of these solutions is the non-analyticity in coupling constant at J = 0. This confirms the non-perturbative nature of the nondiagonal solutions, discussed in the section 2.2.
2. Our numerical study shows that for every exact q = 2 nondiagonal solution there is a numerical real-valued nondiagonal solution in the q = 4 model. It appears that the lack of the zero replicas limit is also present in the interacting case. The solutions are plotted on the Fig.3. Note that in the IR region the diagonal and nondiagonal part exhibit similar behavior up to a numerical factor, and the absolute value of that factor approaches to 1 as M decreases. We also argue that these solutions probably can be analytically constructed in terms of the 1 M expansion, which we will investigate in the future work.

3.
We have shown that the solutions, that we have constructed in the q = 2 case, are suppressed by the diagonal saddle for M > 1 in the path integral (however, they can dominate in the replica partition function, and realize the replica-nondiagonal phase in the M < 1 case). In the q = 4 SYK the replica-nondiagonal saddles are also subleading. This is in agreement with the statement that the SYK partition function is self-averaging at large N [12,16,18,25,38], as verified by exact diagonalization numerics.
We also have studied the analytic nondiagonal solutions in the strong coupling limit βJ ≫ 1. We focused on the class of solutions with the time dependence given by the conformal propagator, and the replica dependence given by a Parisi matrix. The key findings are the following: 4. We were able to obtain valid solutions in the M → 0 limit, using the Parisi ansatz analytic continuation. We have obtained the replica-symmetric solutions, and the solutions with one-step replica symmetry breaking. Meanwhile, the findings of the previous part suggest that these solutions have no well-defined UV completion at M = 0, and therefore have no physical influence on the SYK thermodynamics. Nevertheless, we calculated the regularized free energy on the solutions. Among them there are some that have global minimum regularized free energy value smaller than the replica-diagonal free energy.

5.
We have illustrated on an M = 2 example that the extended reparametrization symmetry in the IR limit has interesting consequences, since it can be used to obtain solutions that have spin-glass-like dynamics for a finite amount of time.
Another solution that can be obtained using reparametrizations is related by analytic continuation to the purification of the replica-diagonal SYK by the thermofield double.
An important technical question that we left unattended in this work is the problem of stability of the replica-nondiagonal solutions that we have obtained in the interacting q = 4 model. The numerical solutions of the exact saddle point equations appeared to be stable to small enough fluctuations in the initial condition. As far as the factorized solutions in the strong coupling regime are concerned, the following comments can be made. We can say that the replica-diagonal solution is stable with respect to the replicadiagonal fluctuations of the field G(τ 1 , τ 2 ), because the quadratic part of the action is determined by the ladder kernel and it has non-negative eigenvalues [3]. The numerical investigation of [20] confirms that there is no instability to general fluctuations of the replica-diagonal saddle point. Since the time dependence of the factorized solutions is the same, we can expect that they also would be stable with respect to the replicadiagonal fluctuations. However, we also can expect that some of the solutions discussed in the section 6 can be unstable to general fluctuations, because they are defined by pairs of complex saddles, and one would need to take into account fluctuations around both saddles for every solution. It could also be instructive to consider other solutions with several steps or continuous replica symmetry breaking to verify these observations. We leave these questions for the future study.
The results of our investigation generally suggest that the common use of the replica-diagonal approximation for study of the quenched averaged quantities in the SYK model is indeed valid at large N , in that we have not found any replica-nondiagonal saddles, contributing to the free energy after the zero replicas limit in the full SYK model. However, the solutions that we have constructed can be interesting in the different aspects.
(i) We have shown that the annealed quantites, which would require a finite number of replicas, do have nontrivial repica-nondiagonal saddle points. In the case of M = 2 and complex-valued β it is shown that non-trivial saddles are crucial for the quantum dynamics of black holes in the work [16,25]. Therefore, an interesting problem is to study the counterparts of our saddles in the spectral form-factor, and to see if they can be responsible for the long-time behavior of this quantity. Definition. The Parisi matrix Q associated with the given tree T with a set D T is defined as follows Q a,a = q 0 , Then we divide m l elements on r l−1 groups with m l−1 elements in each, i.e.
and so on. On the last step we left with r 1 elementary elements, i.e dim r 1 = 1.

A.2 Representation in terms of block matrices I m i
It is convenient to represent the Parisi matrix (A.9) using the family of the block matrices I m i composed on 1's and 0's From this statement follow the Lemma. The space of Parisi matrices with fixed tree T{r 1 , r 2 , r 3 ...r l } is an algebra under both regular and Hadamard matrix products.
Proof. The closeness under the Hadamard product is obvious. For the direct matrix product this follows from representation (A.12) and the following properties (A.17) of I m i matrices corresponding to the same tree. Indeed we define the two Parisi matrices of the rank l (a l+1 = 0): We want to calculate their product, which we denote as We proceed as follows: The four terms come from the multiplication of two brackets between each other. Let us evaluate carefully each term using the relation (A.17). The first term gives The second term yields We relabeled some summation indices and canceled two terms here. The third term can be obtained from T 2 by making the replacement a ↔ b: where we again relabeled the indices in the last term. Therefore, we have presented the four terms as a linear combination of the I-matrices, more specifically four terms = j (U j I m j+1 + V j I m j ) . (A.27) As it will become clear below, this is already proves the lemma. Let us write down the U and V explicitly: As we see, there is a difference between the two: and therefore in R there is a rogue term j D j I m j+1 (A.31) However, we can deal with it using a trivial formula: Thus, we can write the term as follows: Therefore, the lemma is proved. We obtain that W is a Parisi matrix W = j w j (I m j+1 − I m j ) + w 0 I 1 , (A. 34) where its parameters are given by These formulae can be used to directly solve the equation (5.17), if one takes b j = a q−1 j .

A.4 Determinant of the Parisi matrix
The eigenvalues of the Parisi matrix Q are given by the formulae [39] λ 0 = q 0 − q 1 (A.37) This is the result which we use We see that the singular terms at M → 0 limit killed off each other. Taking the limit M → 0, we arrive at the expression log det Q = log(q 0 − q 1 ) + q 1 q 0 − q 1 . (A.42) -49 - The cases of higher l can be calculated analogously. The main idea is to extract the part from the lowest degree bracket which is the same as the next degree bracket, which will always be of the odd degree and carry a singular contribution. Replica-symmetric ansatz. In this case we have two dynamical variables G 0,1 (and corresponding Σ 0,1 (τ ) = J 2 G q−1 0,1 (τ )). To evauate the Pfaffian term in the frequency space we use the formula for the determinant of a Parisi matrix (A.41) for every frequency: . (B.5)