Does the SYK model have a spin glass phase?

We argue that the Sachdev-Ye-Kitaev model has no spin glass phase, based on calculations involving both the nearly-conformal limit and the strongly-coupled Schwarzian limit of the model. This conclusion is supported by numerical computations of eigenvalue statistics with up to 46 Majorana fermions. In addition, we find numerically that the distribution of the ground state energy is Gaussian.


Introduction and summary
The Sachdev-Ye-Kitaev (SYK) model is a disordered quantum mechanical model of N Majorana fermions that has remarkable properties, involving both its direct quantum mechanical description and its holographic dual [1][2][3][4][5][6]. In the large N , strongly coupled limit the model becomes solvable, yet it remains chaotic. It has a master field reformulation that is evocative of a simple bulk description.
While our understanding of its holographic dual is incomplete, many of the model's low temperature properties are reproduced by Jackiw-Teiltelboim gravity [5][6][7][8][9][10]. In particular, the infinite N model has non-zero entropy at zero temperature and a maximal Lyapunov exponent [3,4,11], two properties that are consistent with a bulk description involving an extremal black hole. The SYK model has also found condensed matter applications in strongly-coupled transport and entanglement dynamics [12][13][14][15][16][17][18]. Finally, the fact that the model has a finite-dimensional Hilbert space at finite N allows for straightforward and precise numerical computations.
It is natural to ask whether the model has a transition to a spin glass phase at low temperaturea common occurrence in disordered systems. Indeed, the original Sachdev-Ye (SY) model [1] shares many properties with the SYK model, but in some versions of that model a spin glass transition occurs at a relatively high temperature [19]. A spin glass transition in the SYK model would imply a breakdown of the dual black hole picture. 1 We study this question by looking for two distinct signatures of a spin glass phase. In both cases we find no indication of a spin glass transition, suggesting that the SYK model remains in its well-known paramagnetic phase down to arbitrarily low temperatures.
The first diagnostic is the condensation of replica off-diagonal modes. In the ordinary hightemperature phase, the path integral is dominated by a saddle point that is both diagonal and symmetric in replica space. A deviation from this would indicate a spin glass phase transition. We compute the effective potential for some replica off-diagonal modes in the nearly-conformal limit of the theory, 1 βJ N (here β is the inverse temperature and βJ is the effective coupling; see Appendix A for our conventions). The authors of [19] carried out a similar calculation in the SY model, and for their 'slave fermion' model they found a critical temperature T c Je −c √ N , where c is an order one constant. For the SYK model, a similar comment was made in [5]. We reproduce this estimate of T c using the conformal limit of the SYK model and extend it to arbitrary values of q, the order of the fermion interaction.
Such an exponentially low temperature lies outside the regime of validity of the conformal calculation. Instead, the critical temperature falls within the strongly-coupled Schwarzian limit of the theory, namely 1 N βJ. We repeat the calculation in the Schwarzian theory and find that the effect disappears: The off-diagonal modes we consider are always stable, indicating that there is no spin glass transition. Our analytic results are presented for general q. However, as explained in the text, they are non-trivial only for q ≡ 4 0. The reason is that for the other q values the off-diagonal operator we are considering can never condense.
As a second diagnostic for a spin glass transition, we look for a deviation from Random Matrix Theory (RMT) predictions for the level-spacing statistics at low energies [24]. When the system is in a spin glass phase it loses ergodicity. As a result, we expect its accessible energy states to become uncorrelated, and the level-spacing statistics to no longer follow RMT predictions. In this work we present numerical results for the SYK model with up to N = 46 Majorana fermions and with q = 4. These results were obtained by computing the lowest lying eigenvalues of the Hamiltonian on a cluster of GPUs. Our results are all consistent with RMT predictions, and rule out a spin glass phase for all values of N we tested.
The paper is organized as follows. In Section 2 we carry out the calculation involving the replica off-diagonal modes. In Section 3 we present numerical results for the SYK model, testing RMT predictions involving level-spacing statistics, as well as the distribution of the ground state energy.
Several appendices expand on key points. Appendix A includes our conventions and a brief review of the SYK model. Appendix B includes details of the analytic calculation, and Appendix C describes the numerical methods used in this work. Finally, Appendix D reviews the relation between levelspacing statistics and a spin glass phase in the quantum Sherrington-Kirkpatrick model.

Analytic results
In this section we present an analytic argument against a low-temperature spin glass phase in the SYK model. See Appendix A for a brief review of the model.
In the n-replica theory, the condensation of a replica off-diagonal mode signals a spin glass transition. Such condensation happens when the effective potential of the mode becomes unstable. The effective potential can be computed in the high temperature phase, which is the usual paramagnetic phase described by a replica-diagonal and replica-symmetric saddle point.
In the nearly-conformal limit, we find a predicted spin glass transition (for q = 4) at a temperature T c Je −c √ N with some c > 0. Similar calculations were performed in [25] for the quantum Sherrington-Kirkpatrick (SK) model, and in [19] for the Sachdev-Ye (SY) model. The predicted transition occurs at a temperature βJ N , which is outside the regime of validity of the conformal approximation. Instead, this temperature falls within the strongly-coupled regime of the Schwarzian theory. We repeat the calculation in the Schwarzian theory, and find that the instability actually does not occur.
While these results provide evidence that a spin glass transition does not occur, they do not prove it conclusively. For example, the presence of diagonal, replica-symmetry-breaking solutions may also signal such a transition, and we do not rule out such solutions analytically.

Replica off-diagonal modes
We now introduce the replica off-diagonal modes that will be the focus of the rest of this section, and write down their effective potential to second order in the fields. Let us introduce n replicas, labeled by a, b = 1, . . . , n, and write down the partition function of the replicated theory. After taking the disorder average, we find Let us introduce the Hubbard-Stratonovich field F ab (τ 1 , τ 2 ).
Note that this presentation of the theory in terms of the field F is different than the common presentation in terms of Hubbard-Stratonovich fields G and Σ [4]. The saddle point equation for F is The correlator on the right-hand side is computed in the ordinary SYK theory. We start in the usual (high temperature) phase, dominated by the known replica-diagonal saddle point. In order to detect the putative spin glass phase we will lower the temperature, and look for an instability in modes F ab (τ 1 , τ 2 ) with a = b. If any of these modes condense, that is a signal of a spin glass transition.
For the q = 4 theory, F ab is a 4-fermion operator. Taking a = b, this is the minimal replica off-diagonal operator which can condense. To see why, imagine computing ψ a ψ b where ψ a ψ b is (schematically) a replica off-diagonal operator. Suppose we do this by first computing the fermion path integral, followed by the disorder average. In the first step the replicas are decoupled, and the calculation factorizes as ψ a ψ ψ b ψ J into fermion 1-point functions, which vanish. Here · ψ denotes a fermion path integral and · J denotes disorder averaging.
For the theories with q = 2, 6, 10, . . . , the off-digonal operators F ab involve an odd number of fermions in each replica. The same argument then shows that these operators cannot condense, and we expect their effective potentials to always be stable. For theories with q = 8, 12, . . . an instability in F ab is possible, but our F ab is not the minimal operator that can condense as it involves more than 4 fermions. Therefore, while we carry out the calculation for general q, the resulting evidence against a spin glass transition only applies to theories with q ≡ 4 0 and is strongest for q = 4.
We focus on the time-independent modes F ab with a = b. In Appendix B we compute the quadratic piece in the effective potential of these modes, , and find the following squared-mass.
Here G(τ 1 , τ 2 ) is the usual fermion bilinear operator, defined in Appendix A. The first term on the right-hand side is leading at large N and fixed temperature. A phase transition will happen if m 2 ab becomes negative at sufficiently low temperature. Again, the correlator appearing in (4) is a correlator in the ordinary SYK theory.
Notice that an instability is only possible when q ≡ 4 0, consistent with the argument above. We will assume this from now on. Let us now compute the effective mass in two different limits of the theory.

Nearly-conformal limit
Let us compute the squared mass (4) in the nearly-conformal limit. At leading order in large N the correlator factorizes as G q/2 = G q/2 + · · · . In this limit, the correlator is given by Using this result, we find a log divergence in the integral that we regularize by introducing a cutoff on Euclidean time.
We place the cutoff at ≈ J, where we expect the correlator to become free, G(τ ) free = 1 2 sign(τ ). With these approximations, we find The squared-mass becomes negative at the critical temperature For q > 4, this predicted critical temperature is parametrically smaller than e −N , the typical level spacing. Only the ground state is accessible at this temperature and so this result suggests there is no spin glass phase transition for these theories. For q = 4, we get the predicted transition for q = 4.
Notice that J/T c N , and therefore this critical temperature lies outside the regime of validity of the conformal approximation. The correct description at this temperature is the strongly-coupled Schwarzian theory. We now turn to computing the critical temperature in this limit.

Strongly-coupled Schwarzian limit
The Schwarzian theory is a solvable theory with an inverse coupling constant C = N α S (q)/J [4,[26][27][28][29][30]; see Appendix A for a brief review. In this section we assume that we are in the strong coupling limit C β. We use the results of [29,31] to compute the correlator in the Schwarzian Here c 1 is a constant whose precise value will not be important to us. The factor (β/C) 3/2 exp −2π 2 C/β comes from the normalization by 1/Z [26,31]. We will analyze this formula in two limits.
Let us first consider the regime τ C. The factor exp −(β − τ )k 2 2 /2C is only significant for k 2 C/β 1. This means that to leading order we can set k 2 = 0 in the cosh(2πk 2 ) in the denominator. The k 1 integral is dominated by the range k 1 1, where we can approximate sinh(2πk 1 )/(1 + cosh(2πk 1 )) ≈ 1. Ignoring overall numerical coefficients, we get Note that we recovered the conformal answer.
Next, consider the regime C τ ≤ β/2. Now both τ, β − τ C and the integral in (10) is Let us analyze how the effective mass (4) changes compared to the conformal answer (7). The negative contribution to the effective mass (4) is proportional to 2β Here we used the fact that G(τ ) is symmetric about τ = β/2. Let us compute the τ integral by splitting it into three regions: τ ∈ (0, 1/J), (1/J, C) and (C, β/2). In the first region, the correlator is approximately equal to the free correlator which is a constant, and we get a contribution proportional to 1/J. From the second region we get, using (12), The contribution of the third region is, using (13), We see that the dominant contribution to the τ integral at large N is from the second region, equation (14). Thus from (4), now setting q = 4 for simplicity, we have Here c 2 is a positive constant. Comparing to the conformal answer (7), we see that log(βJ) got replaced by log N . Thus the effective mass is always positive and the mode is stable. The same conclusion holds for other values of q.

Numerical results
In the previous section we studied the condensation of replica off-diagonal modes, which serve as a signature for a spin glass phase transition. For quantum systems, the level spacing statistics are another such signature. Indeed, in an ordinary chaotic system the level spacing statistics obey Random Matrix Theory (RMT) predictions, implying for example level repulsion [24]; in a spin glass the levels are decorrelated and there is no level repulsion. These relations are reviewed in Appendix D for the quantum Sherrington-Kirkpatrick model.
In this section we present numerical results for the spectrum and level spacing statistics of the SYK model with 4-fermion interactions. These results were computed by partially diagonalizing the Hamiltonian, obtaining the energy levels at the edge of the spectrum. Details about the numerical methods used here can be found in Appendix C. Our level spacing results exhibit RMT behavior down to the lowest observed energies, and these results favor our conclusion that the model has no spin glass phase transition at low temperature. In addition, we find numerically that the ground state energy follows a Gaussian distribution. Figure 1 shows the spectral density at the edge of the spectrum. 3 At large N and low energies, the analytic prediction [4,26] is that the density should behave as ρ(E) ∼ √ E − E 0 near the edge. If we simply plot the energy density, we find that there are large fluctuations that mask this effect.

The edge of the spectrum
However, if we shift the energies of each realization by its respective ground state energy (such that In [32], the spectral form factor was introduced as a diagnostic of the late-time dynamics, with connections both to the information paradox and to RMT; see also the recent paper [33]. Our numerical results allow us to test the RMT predictions at larger values of N and at lower temperatures. Figure 2 shows these results. The three notable features discussed in [32], the early 'slope' followed by the late time 'ramp' and 'plateau', are clearly visible. In particular, the ramp is consistent with RMT predictions and indicates a chaotic spectrum.

Level spacing statistics
In order to determine the phase of the system at low energies, we compute the level spacing statistics near the edge of the spectrum. Statistics that agree with Random Matrix Theory predictions imply a chaotic phase, while statistics that follow an exponential distribution (corresponding to uncorrelated energy levels) are a signature of a spin glass phase. See Appendix D for further discussion.
The standard method of computing level spacing statistics involves first 'unfolding' the energy levels such that the mean energy density is one (see for example [24]). This procedure works well in the bulk of the spectrum, but becomes unreliable near the edge due to large fluctuations (such as the ones described in Section 3.1). We compute the level spacing statistics in two different ways that Figure 2: The Spectral Form Factor (denoted g in [32]) with N = 42 Majorana fermions and β = 50, using the same data as in Figure 1.
sidestep this problem. First, we compute the level spacing distribution for the two lowest energy levels, collecting statistics only over different realizations. Second, we compute the distribution of log(r n ) where r n = (E n − E n+1 )/(E n+1 − E n+2 ) for a fixed number of lowest energy states. Both of these distributions can be compared directly with RMT predictions without unfolding [34].
The results, shown in Figure 3, are all consistent with RMT predictions. We compare the computed level-spacing distribution against the Wigner surmise. At the edge of the random matrix spectrum, the eigenvalue density correlations are described by the Airy kernel [35][36][37], while in the bulk of the spectrum they are described by the sine kernel. Despite this difference, it is easy to check empirically that the RMT nearest-neighbor level statistics are well approximated by the Wigner surmise in both cases.
Our results rule out a spin glass phase for the SYK model with up to N = 46 Majorana fermions. The ordinary paramagnetic phase persists down to arbitrarily low energies, even when the thermodynamic approximation breaks down and it is not useful to discuss temperatures.

Ground state energy distribution
In this section we compute the ground state energy distribution of the model (this was previously studied in [38]). The extremal eigenvalues of matrices in common Random Matrix Theory ensembles follow a Tracy-Widom distribution [39,40]. In light of the detailed agreements between SYK and RMT described above, it is natural to ask whether the ground state energy distribution is also consistent with RMT predictions. We observe numerically that this is not the case, and instead   Table 1 lists these results, which show that the Gaussian distribution is clearly preferred.
Next, Figure 5 shows the dependence of the Gaussian parameters on N . We find that the leading large N term in the mean ground state energy is within 10% of the analytic large N prediction [4], which is E 0 ≈ −0.0406N . For the variance we find that a power law provides a good fit.  Table 1: Higher moments for the ground state energy distribution. SYK data was collected from 10 4 realizations for each value of N . GSE data is not shown because the ground state is in the odd charge sector, and we only computed the even charge sector.

A The SYK model
In this section we briefly review some basic properties of the SYK model, following [4]. The dynamical degrees of freedom of the model are N Majorana fermions ψ 1 , . . . , ψ N . The Hamiltonian Here q is a positive even integer, and we introduced the notation ψ i 1 ···iq = ψ i 1 ψ i 2 · · · ψ iq . For each choice of i 1 < i 2 < · · · < i q , the coupling J i 1 ...iq is an independent Gaussian random variable with zero mean and with variance given by The Euclidean time action is The fermion bilinear operator is defined as The nearly-conformal limit of the theory is 1 βJ N . In this limit, the fermion 2-point function is given by This solution defines a replica-diagonal saddle point of the theory, written in terms of its master fields G and Σ [4].
The SYK action has an emergent time reparametrization symmetry τ → f (τ ), which is spontaneously broken by the above solution (21). Furthermore, this symmetry is explicitly broken by corrections to the conformal limit. The effective low-energy action of the theory, which governs the dynamics of the pseudo Nambu-Goldstone modes f (τ ), is the Schwarzian action where C = N α S (q)/J. Here α S (q) is a numerical coefficient whose precise values can be found in [4]. In the weak coupling limit (corresponding to βJ N ), the fluctuations about the saddle point f (τ ) = τ are small, and one can reproduce many results of the SYK model in the conformal limit. In the strong coupling limit (corresponding to βJ N ) the theory is still solvable [26,27,29].

B Derivation of the effective mass
In this appendix we derive equation (4) for the effective mass of the replica off-diagonal modes F ab .
The effective action for these modes is defined by where the replicated partition function was given in (2). Expanding the effective action to quadratic order, we have We show below that only terms in which (a, b) = (c, d) have non-trivial masses, and compute the effective squared-mass m 2 ab ≡ m 2 ab,ab of the time-independent modes. Expanding equations (24) to second order in F ab (τ 1 , τ 2 ), and using (2), we get The second term on the right-hand side can be written as If (a, b) = (c, d) then for some replica (say a) the fermions appear in the correlator all with the same time, and so this correlator vanishes (at leading order) on the replica-symmetric saddle. Therefore only terms where F ab and F cd have the same replicas survive. Let us set a = c, b = d, and consider a specific choice of a, b with a = b.
Let us now compute the correlator appearing in (28) at leading order in large N .
In the last step we kept only the diagonal terms, which give the leading contribution at large N .
Plugging this in (28) and then (26), we find Focusing on the time-independent mode of F ab , we get equation (4) as advertised.

C Numerical methods
In this appendix we provide details about the numerical methods used to compute the results of Section 3. We used two independent implementations to test our results, one running on CPUs and one on GPUs. The GPU implementation can be found at https://github.com/guygurari/syk.
with properly chosen O a L/R operators and g a couplings. Computing the tensor product operations in the above representation will bring us back to the standard approach to exact diagnoalization.
However, tensor product operations are quite expensive computationally and storing the resulting huge matrices is costly. It is possible to avoid doing such unnecessary costly operations and still perform diagonalization algorithms efficiently.
The Lanczos algorithm is one of the most popular methods for obtaining low energy eigenvectors and eigenvalues of sparse matrices [41]. It is a power iteration method based on successive matrixvector multiplication operations, v b+1 = Hv b , starting from an initial random vector v 0 . The resulting v b 's form basis vectors for the Lanczos diagonalization procedure. The desired vector operations can be implemented more efficiently using the Schmidt decomposition of v b vectors, where a,L/R are orthogonal basis vectors defined on the L (R) subsystem. We then utilize the following unitary (duality) transformation on the right side : v Next, we consider v b+1 = Hv b . It can be verified that v b+1 = reshape (v b , D L , D R ) can be evaluated This way we never need to explicitly compute the tensor products O a L ⊗ O a R , and instead we just need to store the when O a L and O a R operators are dense matrices. For sparse operators, the time complexity is unaffected by the above scheme. However, we have noticed that in practice the complexity can drop significantly using the above method, especially for long range Hamiltnonians such as SYK (indeed, the overall prefactor of the time complexity decreases).
We implemented the Lanczos algorithm as described on GPUs, taking advantage of their ability to carry out highly parallel calculations. In common implementations of the Lanczos algorithm one keeps track of previously computed eigenvectors, in order to overcome the inherent numerical instabilities of the algorithm. Despite the lower space complexity described above, this method is still too costly to run on GPUs due to their relatively limited RAM. Instead, in the GPU code we used an alternative implementation of the Lanczos algorithm which does not need to keep track of the eigenvectors [41]. This method is useful when one is only interested in the eigenvalues of the matrix.

D The Quantum Sherrington-Kirkpatrick model
In this appendix, we present numerical calculations of the eigenvalue statistics in the quantum Sherrington-Kirkpatrick model [42]. The Hamiltonian is that of the transverse field Ising model on N sites, but with random infinite-range couplings: where X i and Z i are the Pauli-X and Pauli-Z matrices on site i. The sum in the first term runs over all pairs in the system, and the couplings J ij are independent Gaussian random variables with mean zero and variance J 2 ij = 1/N . This Hamiltonian has a Z 2 symmetry represented by the unitary operator U = Z 1 . . . Z N .
When Γ = 0 in the Hamiltonian (35), all the terms in the Hamiltonian commute and the model reduces to the classical Sherrington-Kirkpatrick model, which is well-known to have a spin glass phase [43][44][45] at low temperatures. The spin glass phase persists at small Γ. The spin glass phase can be destroyed by either increasing Γ or by increasing the temperature beyond their critical values.
A cartoon phase diagram of the model is shown in Figure 6.

Spin Glass
Paramagnet Transverse field Γ Temperature T Figure 6: Schematic phase diagram of the quantum Sherrington-Kirkpatrick model. There is a spin glass phase at small external field Γ and small temperatures.
We project the Hamiltonian to the even Z 2 sector and perform exact diagonalization on a system of N = 12 spins. The level spacing statistics calculation was described in Section 3.2 (see also [34]).
In Figure 7, we show the distribution of log r n when Γ = 0.1. At this value of Γ, the low temperature phase is a spin glass. Thus the low energy part of the spectrum should exhibit exponential statistics, as is clearly visible in the left panel of Figure 7. We use 200 disorder realizations and the lowest 50 states from each realization. The high temperature phase is ergodic, and thus states drawn from the middle of the spectrum should exhibit GOE statistics. This is also clearly visible in the right panel of Figure 7. Here we take 200 disorder realizations and the middle 50 states from the spectrum of each realization.
Finally, consider the case where Γ is large. Here there is no spin glass phase at any temperature, so even the low energy part of spectrum should exhibit GOE statistics. This is confirmed in Figure 8.  Figure 7: Distribution of log r n when Γ = 0.1 for low lying states (left) and for states in the middle of the spectrum (right). The blue dots are numerical data, the red curve is the exponential distribution, and the green curve is the GOE ensemble prediction. There is a spin glass phase at low temperatures, and consequently the distribution is exponential. Since the spin glass phase is wiped out at high temperatures, the states from the middle of the spectrum follow GOE statistics. There is no spin glass phase and thus even the low-lying states follow GOE statistics. The blue dots are numerical data, the red curve is for the exponential distribution, and the dark green curve is for the GOE ensemble.