Quenched vs Annealed: Glassiness from SK to SYK

We show that any SYK-like model with finite-body interactions among \textit{local} degrees of freedom, e.g., bosons or spins, has a fundamental difference from the standard fermionic model: the former fails to be described by an annealed free energy at low temperature. In this respect, such models more closely resemble spin glasses. We demonstrate this by two means: first, a general theorem proving that the annealed free energy is divergent at low temperature in any model with a tensor product Hilbert space; and second, a replica treatment of two prominent examples which exhibit phase transitions from an"annealed"phase to a"non-annealed"phase as a function of temperature. We further show that this effect appears only at $O(N)$'th order in a $1/N$ expansion, even though lower-order terms misleadingly seem to converge. Our results prove that the non-bosonic nature of the particles in SYK is an essential ingredient for its physics, highlight connections between local models and spin glasses, and raise important questions as to the role of fermions and/or glassiness in holography.

What the SYK model does not exhibit is spin glass physics [25][26][27]. This is surprising, because the SYK Hamiltonian bears a striking similarity to the quintessential mean-field models of spin glass theory [28][29][30]. Both the SYK and spin glass models are defined by randomstrength interactions among all degrees of freedom. Here we show that the essential difference is the fermionic nature of the particles in SYK: any model with strictly local degrees of freedom will share much more in common with spin glasses.
This result is relevant because interest in SYK physics has spread to generalizations of the original model. To name a few: including multiple flavors of fermions [31], using bosonic particles [2,32,33], using spins [34,35], forming lattices of SYK models [36,37], and introducing supersymmetry [38,39]. With the analysis presented in this paper, we are able to immediately identify large classes of such models in which the potential for glassiness must be carefully addressed.
On the spin glass side, all-to-all disordered models have featured prominently for decades. Sherrington and Kirkpatrick first introduced a system of Ising spins with infinite-range random interactions which exhibits an intricate spin glass phase [28,40,41]. The model has been extended in numerous directions, both classical and quantum, many of which are central to the field in their own right: p-body interactions [29,42,43], spherical spins [44][45][46], Potts spins [47,48], Heisenberg interactions [1,2,30,49], and transverse fields [50][51][52], among many others. These variants all share certain phenomena which unite them as spin glasses. As one lowers the temperature, the system first undergoes a "dynamical" transition at temperature T d , below which dynamical correlation functions never fully decay. The system experiences a further "static" transition at a potentially lower T s , below which one can detect frozen magnetization patterns in the equilibrium Gibbs distribution. Certain systems undergo a third "Gardner" transition at an even lower T g , below which the magnetization patterns become more complex, with sub-patterns and so on. For pedagogical expositions of the physics, see Refs. [53][54][55].
All indications are that the SYK model does not show any such behavior [1,[25][26][27]56]. This raises multiple questions, chief among which is simply: which generalizations of the SYK model do have spin glass phases? Presumably such glassiness would rule out any connection to quantum gravity (although that is itself an important open question). It has long been known that the bosonic variant of SYK is a spin glass [2,3,32], yet this is merely one model out of the multitude which could arise. A recent numerical study on small systems found evidence suggesting that the hard-core bosonic variant is a spin glass as well [33]. Beyond this, the question has remained unexplored. There has been no general framework for understanding when all-to-all disordered systems behave as spin glasses rather than SYK.
This paper aims to fill that gap. On a technical level, generic models can be analyzed in two ways, and we address both. The first relies on the replica formalism: one expresses the moments of the partition function as a path integral and uses standard mean-field techniques to obtain the free energy [53][54][55]. One can circumvent replicas by making the "annealed" approximation, namely replacing the partition function by its first moment at the outset. The second approach is to organize the diagrammatic expansion of the propagator in powers of system size N . One averages each term over the disorder and finds that a summable set of diagrams (the so-called "melons" in SYK) gives the leading-in-N contribution. This ultimately gives the same results as the annealed approximation. Even though the annealed approximation appears to be correct for the SYK model, it is in general extremely unreliable at low temperature. Indeed, breakdown of the annealed approximation is often what signals entry into a spin glass phase. We will be studying this breakdown and its consequences in generic all-to-all disordered systems.
In Sec. II, we introduce our notation and the specific models which will serve as our examples. In Sec. III, we prove that the annealed approximation cannot hold at low temperature in any model for which the Hilbert space is a tensor product. This includes bosons (soft-& hardcore), spins, distinguishable particles, etc. It shows that all such models are fundamentally different from SYK. In Sec. IV, we then give a more detailed and transparent analysis of the hard-core bosonic and quantum p-spin models. Despite the models not being fully solvable, we show that each undergoes a transition from an annealed phase at high temperature to a non-annealed phase at low temperature. Lastly, in Sec. V, we use a concrete example to demonstrate the difficulty in obtaining such results through a 1/N expansion.

II. MODELS AND DEFINITIONS
Here we define the models of interest, starting with the original SYK model and then introducing various modifications. All of the models discussed here are in fact ensembles of Hamiltonians given by Gaussian random couplings. We also give a very brief description of the replica method. More detailed accounts can be found in the references.
• Fermionic models: The original SYK model is defined using N Majorana (i.e., Hermitian) fermion operatorsγ i . Note that the Hilbert space of the theory has dimension 2 N/2 . The Hamiltonian, which has an even integer q as a parameter, is where the couplings J i1···iq are independent Gaussian random variables with mean zero and variance One can also consider the analogous complex SYK model, where the Majorana operators are replaced by complex fermionsĉ i andĉ † i (p ≡ q/2 of each): Here and throughout, we use a convenient notation in which the multi-index I represents a set of p indices i 1 < · · · < i p arranged in increasing order.
Thus H cSYK consists of all possible p-body interactions. The couplings J II are again independent Gaussians, but now complex with and such that J I I = J * II . One reason for considering H cSYK as opposed to H SYK (or vice-versa) is that H cSYK has a conserved particle number.
• Bosonic models: The bosonic SYK model simply replaces the fermionic operatorsĉ i with bosonic operatorsb i : The couplings J II remain exactly as in Eq. (4). An issue with this definition is that in the grandcanonical ensemble, where the number of particles is unlimited, H bSYK is unbounded from below. One could therefore work at fixed particle number, as past works on the bosonic SYK model have done [1,2,32], or one could interpret theb i as hard-core bosons [33] (i.e., exclude double occupancies on sites). Either choice guarantees that the model has a definite ground state. We shall do the latter: in addition to being more interesting (in the sense that much less is known about it), the hardcore model has the benefit of having a well-defined grand-canonical ensemble.
• Spin models: The quantum p-spin model consists of all-to-all p-body interactions among spinŝ σ α i (α ∈ {x, y, z}): where I is the same multi-index as before and A = {α 1 , · · · , α p }. We will use spin-1/2, but our results apply to any spin. The couplings J A I are real Gaussians with variance Connections between this spin model and SYK have recently been explored in Refs. [34,35].
It should be stressed that our conclusions are in no way restricted to these models. We focus on those listed here solely for the sake of concreteness and current relevance. Regardless of the model, one is always faced with the question of how to treat the random couplings (the "disorder"). Assuming the ultimate goal is to calculate the statistics of physical observables, an important quantity is the "quenched" free energy: where E[ · ] denotes the average over random couplings and Tr[ · ] is the usual sum over states. Derivatives of f (β) clearly give the disorder-averaged values of observables, exactly as the free energy does in non-random systems. f (β) is extremely difficult to evaluate, even for classical systems. The replica method is one of the few ways to make analytic progress. It is based on the identity E ln Tre −βH = lim One evaluates the average on the right-hand side for integer n, interpreting (Tre −βH ) n as the partition function for n uncoupled "replicas" of the system, i.e., where {Ψ} is a complete set of states. In the first line, the operators and states live in the original Hilbert space H, whereas in the second line, they live in the product space H ⊗n . Assuming one can obtain an analytic expression for the disorder average of Eq. (10), one then pretends that n is an arbitrary real number and takes the n → 0 limit. This technique is clearly not rigorous. It has nonetheless been tremendously successful in the study of disordered systems [53][54][55].
A drastic but useful approximation which avoids replicas entirely is to interchange the disorder average and logarithm in the definition of the quenched free energy (and then take the average inside the trace). This gives the "annealed" free energy: Note that derivatives of f (ann) (β) do not correspond to physical quantities. One often finds that f ∼ f (ann) at high temperature but that f (ann) gives patently incorrect results at low temperature (see Sec. III). The SYK model seems to be the only known non-trivial counterexample.
As an aside, the terminology "quenched" versus "annealed" comes from metallurgy, and refers to whether fluctuations in the disorder (accounted for by E[ · ]) are treated on the same footing as thermal fluctuations in the degrees of freedom (accounted for by the trace). Eq. (8) treats the disorder as fixed when computing observables and only afterwards averages over disorder, whereas Eq. (11) sums over fluctuations in both simultaneously.

III. BREAKDOWN OF THE ANNEALED APPROXIMATION IN TENSOR PRODUCT MODELS
Here we prove a general result: the annealed free energy cannot be correct at low temperature for any all-toall model with a tensor product structure. Specifically, consider any N -particle Hamiltonian of the form where I denotes sets of p particles and A denotes sets of p indices from some group of size k, and J A I is Gaussian with We shall take the operators O α i to be Hermitian, but models such as complex SYK involving non-Hermitian operators can be treated in the exact same manner. The only restriction we place on the operators is that they obey a tensor product structure: the Hilbert space H is a tensor product H 1 ⊗ · · · ⊗ H N andÔ α i is shorthand for 1 1 ⊗ · · · ⊗Ô α i ⊗ · · · ⊗ 1 N . The quenched and annealed free energies are, respectively, We prove that there is a finite β * such that for β > β * , A. Warm-up Let us first consider a classical model, for which the annealed free energy is easily computed. The Sherrington-Kirkpatrick (SK) model mentioned in the introduction is with Ising spins σ z i and Var[J ij ] = 1/N . A simple calculation gives and thus Yet if Eq. (18) were the correct expression for the average free energy, then the average energy per spin would be = −β/2 and the average entropy would be s( ) = ln 2− 2 . This cannot be, since the entropy in a discrete configuration space is non-negative: the number of configurations Ω( ) within a small energy window around is a non-negative integer, thus lim N →∞ N −1 ln Ω( ) is either −∞ or non-negative. The annealed free energy of the SK model must be invalid for < − √ ln 2, i.e., β > 2 √ ln 2.

B. Generic tensor product models
The statement that the entropy must be non-negative applies equally well to quantum systems: simply replace the word "configurations" by "energy eigenstates". For any Hamiltonian H g of the form in Eq. (12), we give an upper bound to f (ann) which diverges to −∞ as T ≡ 1/β → 0. It follows that the annealed entropy, being −∂f (ann) /∂T , must diverge to −∞ as T → 0, and thus cannot be correct below a certain temperature. See Fig. 1 for a sketch of the situation.
Since the variousÔ α may not commute for different α, we cannot directly evaluate the annealed free energy as for the SK model. Yet we always have Jensen's inequality: for any quantum state |Ψ . Summing Eq. (19) over a complete set of states, averaging over disorder, and taking the logarithm, we find that Note that Eq. (20) holds for any basis |Ψ used on the right-hand side.
The tensor product structure allows us to use a product basis, i.e., Furthermore, since the operatorsÔ α are not identically 0, there must be some single-particle state |ψ * i for which at least for some α. Use this |ψ * i as a basis state. Then where the omitted terms (coming from A = {α, · · · , α}) are positive. Inserting into Eq. (20) gives our final bound: Clearly f (ann) → −∞ as β → ∞, as claimed.
The divergence of f (ann) at low temperature is not limited to Gaussian disorder. The Gaussian coupling distribution was used only to evaluate the average in Eq. (20), and an analogous bound can be obtained for any other distribution. For example, suppose each J A I has some alternate probability density P (J) for which the mean is zero and the variance is still given by Eq. (13). Assume P (J) falls off faster than exponentially for J 2 Var[J], so that we can safely expand inside the average: We can proceed with the proof as before and obtain the same Eq. (24), for any such P (J).
For the sake of concreteness, we next consider some specific models.

C. Example: Quantum p-spin
A natural basis to use for the quantum p-spin Hamiltonian (Eq. (6)) is theσ z i eigenstates | ↑ and | ↓ . Both states have expectation values | σ z (20) gives the bound The extra term compared to Eq. (24) comes from summing over the 2 N basis states, which we neglected for simplicity in the general treatment.

D. Example: Hard-core bosonic SYK
In this case (Eq. (5)), we have that

E. Example: Complex fermionic SYK
Note that the bound in Eq. (20) always applies, regardless of whether the Hilbert space is a tensor product or not. In particular, it holds for the fermionic SYK model (both real and complex), for which the annealed free energy seems to be correct at all temperatures. It is informative to see how this result is consistent with the bounds obtained through Eq. (20), in contrast to the examples above.
Given the similarity between hard-core bosons and fermions, and given that the states (|0 ± |1 )/ √ 2 yielded a free energy diverging at low temperature in the former, let us consider the analogous basis in the fermionic Hilbert space: Here the state index Ψ ∈ {0, 1} N is denoted as a vector ψ, and similarly for s. Starting from Eq. (20), we need to evaluate Note that to leading order in N , none of the i k are equal to any of the j l . A given term vanishes unless s i1 = · · · = s ip = 0, s j1 = · · · = s jp = 1. Furthermore, we need s k = s k except for k ∈ {i 1 , · · · , i p , j 1 , · · · , j p }, in which case s k = 1 − s k . Thus (−1) ψ·( s+ s ) = (−1) ψi 1 +···+ψi p +ψj p +···+ψj 1 , which can be taken outside the sum.
Additional minus signs come from rearranging the fermion operators. First note that factors ofĉ i1ĉ † i1 ,ĉ j3ĉ † j3 , etc., in which the two matching operators are adjacent, can be treated as the identity: as a pair, they commute past all other operators, andĉ iĉ † i |0 = |0 . Thus owing to the initial order of theĉ † l in Eq. (29) (the index increases from left to right), one can convince oneself that we obtain a factor of −1 for each k less than i 1 , each k less than i 2 , each l less than j 1 , each l less than j 2 , and so on. But now suppose that j 1 ≤ i 1 − 2. Each choice of s can be associated with an r according to s j1−1 = 1 − r j1−1 , s j1+1 = 1 − r j1+1 , with all other s l = r l . The two vectors give contributions differing by exactly one minus sign, and therefore sum to 0.
The lesson is that Eq. (30) evaluates to 0 unless every i & j index is adjacent to another, e.g., i 1 = j 1 − 1 or i 2 = j 1 + 1. Yet this restricts the number of free indices for us to sum over, and we needed all 2p to be free in order to obtain an extensive bound on f (ann) (a factor of N 2p to compensate for N −(2p−1) from the coupling variance). In the thermodynamic limit, the only bound we obtain in this case is The right-hand side does not diverge as β → ∞, and f (ann) cSYK has the potential to remain correct even at zero temperature.

IV. REPLICA ANALYSIS FOR SPECIFIC MODELS
Much additional insight comes from considering the replica analysis in detail for specific models. We shall focus on the hard-core bosonic and p-spin models. Yet keep in mind that even though we limit ourselves to these two for ease of presentation, our analysis is in fact much more general. It can be applied with minimal modifications to any model which admits a path integral representation, even if local constraints on the fields are required.
Furthermore, the replica analysis allows us to make conclusions about the high-temperature behavior of the models. Indeed, we show that the hard-core bosonic and p-spin models undergo genuine phase transitions: for each, there exists a β c such that the free energy equals f (ann) for β < β c and does not for β > β c . We are able to say this without needing to calculate the precise functional behavior of f (ann) .

A. Hard-core bosonic SYK
The hard-core bosonic SYK model is given by Eq. (5), reproduced here: To construct a path integral representation of the partition function, we express each hard-core boson operator b i as a pair of fermions together with a constraint: The partition function is then with Note that we enforce the constraint by way of a Lagrange multiplier µ i on each site. Also note that we use a slightly non-standard definition of imaginary time: τ ranges from 0 to 1 for all β. This will be convenient in what follows. As described in Sec. II, we now evaluate E[Z n bSYK ]. To save space, we give the steps of the calculation in Appendix A. The method is standard, and analogous calculations can be found in, e.g., Refs. [3,51,54,55]. The result is a path integral over an order parameter G rr (τ, τ ) and Lagrange multiplier F rr (τ, τ ): where The indices r and r denote different replicas: r, r ∈ {1, · · · , n}.
In the thermodynamic limit, Eq. (36) is dominated by the saddle-point value, whose location is determined by the equations where · eff denotes an expectation value using the effective action of Eq. (38). From the analysis in Appendix A, one can show that in physical terms, the solution G rr (τ, τ ) to Eqs. (39) and (40) is simply the equilibrium Green's function for the hard-core bosons in the original model: where T denotes time ordering.
Thus far, all calculations have been exact. We cannot proceed any further in full generality, since (unlike in the SYK model) the remaining action for h and a is not quadratic. Nonetheless, we shall use this starting point to both confirm that the annealed free energy diverges at low temperature and show that it is correct at high temperature.

Low temperature
If we set n = 1 in Eq. (36), then we in fact have an expression for the annealed free energy: In terms of the order parameter G(τ −τ ) (note the lack of replica indices and that we have assumed time translation invariance), the expression for Φ 1 is and the saddle-point equations are At low temperature, the maximizer of Eq. (43) is static, i.e., independent of τ . We show this self-consistently. Given a τ -independent F , we can perform a Hubbard-Stratonovich transformation on the remaining path integral: The path integral on the second line is precisely that of a single hard-core boson with z-dependent Hamiltonian We at last have a tractable expression. Eq. (45) becomes The right-hand side must be independent of τ in order to be consistent, and this is indeed what happens at low temperature: for τ 1/β ∼ 0, we have Returning to Eq. (43), the annealed free energy is where the ellipses denote terms subleading in β.
Of course, F (τ ) is not strictly static at low but nonzero temperature. Rather, it has a static component F and a correction ∆F (τ ), where the correction is nonnegligible only for τ 1/β. The presence of ∆F (τ ) does not change the fact that the correlation time is O(1/β), and thus this ansatz for F (τ ) is fully self-consistent. Furthermore, ∆F (τ ) only gives subleading corrections to f (ann) . The expression shown in Eq. (50) is correct to leading order.
Finally, note that while we obtained the annealed free energy by setting n = 1 in Φ n [G, F ] (Eq. (37)), the same expression results from instead setting all inter-replica order parameters to 0: set G rr = F rr = 0 for r = r , and take the n → 0 limit as prescribed (see Sec. II). Since we know that this expression cannot be correct at low temperature, it follows that the true equilibrium value of the order parameter, whatever it may be, cannot be diagonal in replica indices.
The conclusion is that in the hard-core bosonic SYK model, the autocorrelation function G(τ ) develops a static component as T → 0 which is responsible for the divergent annealed free energy. This in turn implies that the system must no longer be replica-diagonal. It also suggests why the fermionic model should behave differently: there, G(τ ) cannot have a non-zero static component because the Fourier transform only has weight on odd multiples of π.

High temperature
To show that the annealed free energy is correct at high temperature, we place a bound on the probability of a random disorder realization having free energy other than f (ann) . Write the (random) partition function Z as ZE[Z]. Chebyshev's inequality states that We shall show that for β less than a certain value, it follows from Eq. (51) that f bSYK = f (ann) bSYK with probability approaching 1 in the thermodynamic limit.
The second moment of Z bSYK is obtained by setting n = 2 in Eq. (36). We have two order parameters: the intra-replica correlator G rr (τ, τ ) ≡ G(τ −τ ) (r ∈ {1, 2}) and the inter-replica correlator G 12 (τ, τ ) ≡ Q, which we take to be static. Thus The saddle-point equations are We immediately have one solution to the saddle-point equations. Denote the solution to the annealed equations, Eqs. (44) and (45), by G eq (τ ) and F eq (τ ). It is self-consistent to set in Eqs. (55) through (58). Note that this is not a trivial statement: it relies on the fact that in the absence of any inter-replica coupling, each replica has a separate U (1) symmetry which ensures h r (τ ) * a r (τ ) = 0. If the action were to include an explicit U (1)-breaking term, then Q = 0 would not be a valid solution.
The question is now whether Q = 0 is the dominant solution to the saddle-point equations, i.e., that which maximizes Φ 2 . To address this carefully, in Eq. (54), denote every part of the expression except for Q 2p and λQ by A 2 [G, F ; λ]. This lets us write the second moment of Z bSYK as follows: This is easiest to see starting from Eq. (A4) in Appendix A. As a result, we can write We have already established that Λ 2 [0] = 0, by virtue of Q = 0 satisfying the saddle-point equations. Furthermore, evaluating Eq. (61) by saddle-point demonstrates that Λ 2 [Q] is the Legendre-Fenchel transform of A 2 [G eq , F eq ; λ] and must therefore be convex [57]. Thus Q = 0 is the unique minimum of Λ 2 [Q]. Finally, a direct calculation starting from Eq. (61) gives i.e., the curvature of Λ 2 remains non-zero even as β → 0. These observations imply that β 2 Q 2p −Λ 2 [Q], although not itself concave, has its global maximum at Q = 0 for β less than a certain non-zero value. Furthermore, for such β, This establishes what we claimed: there exists a critical temperature above which Var[Z bSYK ]/E[Z bSYK ] 2 → 0 in the thermodynamic limit and f bSYK = f (ann) bSYK .

B. Quantum p-spin
The quantum p-spin model is given by Eq. (6), reproduced here: Our treatment of it will be very similar to that of the bosonic SYK model, and we include it here to highlight the generality of the method. For that reason, we will present only the results of each step, and leave the details to be filled in by analogy with Sec. IV A. We express the partition function in terms of spin coherent states |Ω . In fact, the only features of the states that we need are the identities [58,59] (67) The integrals are over the unit sphere, and Ω x i = sin θ i cos φ i , etc. The partition function is (68) In this notation, we are including the overlaps between coherent states in the integration measure, i.e., 2π We will never need to express the overlaps in continuum notation.
The replicated, disorder-averaged partition function is The saddle-point equations are then The notation is the same as before: r and r are replica indices, and · eff denotes an expectation value using the effective action of Eq. (72).

Low temperature
The n = 1 action is We again make a static ansatz for F (τ ), which will turn out to be consistent at low temperature. A Hubbard-Stratonovich transformation gives The remaining path integral is that of a single spin-1/2 in a magnetic field proportional to h, thus we can evaluate the Ω(τ ) · Ω(τ ) correlator directly: which is simply G ∼ 1/27 in the limit β → ∞ with τ 1/β. Finally, the annealed free energy is which indeed diverges at low temperature.
High temperature Using the same notation as in Sec. IV A, the n = 2 action is One saddle-point of Φ 2 is at where G eq (τ ) and F eq (τ ) are the order parameters which maximize Φ 1 . Without any coupling between the two replicas in Eq. (79), Ω α 1 (τ )Ω α 2 (τ ) = Ω α 1 (τ ) Ω α 2 (τ ) , and Ω α r (τ ) = 0 owing to the (statistical) symmetry of the original Hamiltonian. By establishing that this is the dominant saddle-point at high temperature, we show that with probability 1, exactly as done for the hard-core bosonic model.
where A 2 consists of all the remaining terms in Φ 2 . By the same arguments as in Sec. IV A, we have that: (83) Thus for β less than some non-zero value, E[Z 2 p ] is dominated by Q = 0 and The free energy then agrees with the annealed value [60].

V. THE DANGER IN 1/N EXPANSIONS
The preceding sections have been in the spirit of the replica formalism, but there is another technique for studying all-to-all disordered models: expanding the free energy in powers of system size N . Many studies of the SYK model and its variants have taken the latter approach [31,35,61]. Although in principle one could obtain all of the above results through a 1/N expansion, the correct low-temperature physics cannot be identified without taking subtle issues of convergence seriously. The purpose of this final section is to present an example of the issues which arise, as an argument in favor of the replica method over 1/N expansions.
First consider the structure of a 1/N expansion, say for the quantum p-spin model for concreteness. Suppose one wishes to compute the moments of the partition function, E[Z n p ]. We can expand the exponentials: (−β) L1+···+Ln L 1 ! · · · L n ! Tr 1 · · · Tr n E H L1 p,1 · · · H Ln p,n .
Note that in the second line, the n replicas are considered as separate degrees of freedom, each with its own trace. However, the L r factors of H p,r all involve the same spins of replica r, and every factor contains the same Gaussian couplings J A I . Since H p is linear in the couplings, the product H L1 p,1 · · · H Ln p,n is a sum of products of Gaussians, and the disorder average is given by all pairwise contractions according to Wick's theorem. These features are all naturally expressed in terms of chord diagrams [34,35], which we describe in Appendix B. Evaluation of Eq. (85) is then reduced to a sum over chord diagrams. Each diagram comes with a power of N , which allows the sum to be organized as a 1/N expansion.
We further show in Appendix B that, assuming all L r N , the diagrams having contractions between replicas are subleading. In other words, the disorder average factors to leading order: Since in the thermodynamic limit L r N for any fixed L r , the naive conclusion would be that the entire sum factors, and thus E[Z n p ] ∼ E[Z p ] n . In particular, The error is in assuming that the dominant terms of Eq. (85) have L ∼ O(1). Since we expect the energy to be extensive, i.e., H p ∼ O(N ), the expansion of e −βHp should be dominated by L ∼ O(N ). We must at the very least include such L in our evaluation of Eq. (85).
The non-commutativity of the operators in the quantum model makes it difficult to be any more quantitative. Thus instead consider the simpler classical model: where the sum is again over all multi-indices of p spins, and Var[J I ] = p!/2N p−1 . Note that p = 2 is precisely the SK model described in Sec. III (Eq. (16)). Every statement made above about the quantum pspin model can also be made about the classical model, and in the classical model we can confirm our suspicion that the breakdown of the annealed approximation appears only at O(N )'th order in the 1/N expansion. We start with A term of Eq. (88) with given (L 1 , L 2 ) contains L 1 + L 2 factors of the Gaussian couplings, which are contracted in pairs. Some contractions will connect spins on the first replica to spins on the second. By organizing the expansion in terms of the number L of such pairings, as detailed in Appendix B, we have that where in the first line σ z Ij r ≡ i∈Ij σ z ir (r is the replica index), and in the second line the nested sum has been factored. We also used that E[Z cl ] = e N (ln 2+β 2 /4) .
To proceed, write the trace as with The integral over Q can be evaluated by saddle point, leaving us with a single sum over L.
First consider L ∼ O(1) with respect to N . The saddle point is at Q = 1/2, and the integral goes as (pL/2 − 1)!!N −pL/2 . We have that Every term in the sum over L is subleading in N , for any β (except if p = 2, in which case see [60]). This would seem to say that E[Z 2 cl ] ∼ E[Z cl ] 2 , even though that cannot possibly be the correct result at low temperature.
However, consider L ∼ O(N ). The saddle point Q * is now given by the equation where l ≡ L/N , and the sum over L can be written (ignoring sub-exponential prefactors)  Yet if β is large, the maximum of g(l) is positive. The fluctuations in Z cl become greater than the mean, and we can no longer claim that the annealed free energy is correct.
We have shown that a 1/N expansion of the partition function (and thus of the free energy) converges at small β but diverges at large β. Furthermore, note that the saddle point of Eq. (94) is at l * ∼ β 2 /2 for large β, i.e., L * ∼ N β 2 /2. Were one to take the N → ∞ limit be-fore resumming the series, one would miss the divergence entirely, and indeed overlook much of what makes these spin glass models interesting.

VI. CONCLUSION
We have demonstrated that the annealed approximation breaks down at low temperature in any all-to-all disordered model with finite-body interactions and a tensor product Hilbert space. This encompasses many in the family of SYK-like models, such as the bosonic variants and the quantum p-spin model. Furthermore, we have shown that, at least in the hard-core bosonic and quantum p-spin models (although the technique can easily be generalized), the partition function is self-averaging at high temperature. Thus we have identified two distinct phases: one in which the free energy equals the annealed value, and one in which it does not. These results were obtained using rigorous bounds on the annealed free energy and the replica technique. Note that we did not rely on any of the more cryptic aspects of the replica method (taking the number of replicas to 0 and maximizing rather than minimizing the free energy). Finally, we have highlighted the subtleties that come with applying 1/N expansions to such models.
Strictly speaking, these results are not enough to prove that the models are spin glasses at low temperature. Spin glass order is characterized by an overlap matrix in which the permutation symmetry is broken ("replica symmetry breaking"), whereas we have shown only that the matrix cannot be diagonal. In more physical terms, a spin glass has multiple low temperature states, whereas we have shown only the existence of some low temperature state distinct from the high temperature state.
That said, the results established here do force one to confront the issue of glassiness. The standard annealed approximation cannot accurately describe the effects of disorder in any tensor product model, and one must use an approach which allows for non-diagonal and potentially symmetry-broken replica order parameters. In particular, this statement applies to many models of current interest in the context of SYK physics. Whether replica symmetry is broken or merely non-diagonal in any specific model is an interesting open question which requires further analysis.
As for the relevance of these models to holography, it is still possible that some might have gravitational duals despite the breakdown of the annealed approximation. The precise dynamics cannot be exactly as in fermionic SYK, since that model is described by the annealed approximation, but a more complex gravitational theory is not ruled out. It is also possible that glassiness and gravitational dynamics can coexist in an interesting way, e.g., Refs. [62][63][64]. These are all important questions that remain to be investigated.
There is one potential way for the annealed free energy to remain accurate at low temperature even in tensor product models: have an interaction degree which increases with system size. Note that every bound obtained here no longer diverges if the p → ∞ limit is taken before the T → 0 limit. This does not prove that the annealed approximation holds, but we cannot claim that it must break down in such models. One example is the double-scaling limit studied in Refs. [34,35,65,66], where p ∼ √ N . It was argued that the quantum p-spin model has a Schwarzian density of states in this limit. In view of our results, it would clearly be desirable to have a more detailed understanding of the low-energy physics for general p.
Finally, it is interesting to note that every system currently known to have a simple gravitational dual includes fermionic degrees of freedom. This could be a streetlight effect, perhaps related to the difficulty of reliably studying non-supersymmetric theories at strong coupling. However, here we have uncovered a general result preventing a wide class of bosonic theories from exhibiting the simplest kind of gravitational dynamics known to occur in a corresponding fermionic theory. Perhaps this is one example of a general class of constraints which places purely bosonic theories of gravity into the swampland [67].
We present the details for the hard-core bosonic model. The quantum p-spin model and others proceed analogously. Beginning from Eq. (33) The sums over spins/multi-indices now come alongside sums over replica indices r, r ∈ {1, · · · , n}. Note that, to leading order in N , · h i p r (τ ) * a i p r (τ )a i p r (τ ) * h i p r (τ ) · · · h i 1 r (τ ) * a i 1 r (τ )a i 1 r (τ ) * h i 1 r (τ ) In the large-N limit, the path integral is dominated by a specific value of G rr (τ, τ ). We determine this saddle point by introducing a Lagrange multipler F rr (τ, τ ), and thus have where S (eff) is given by Eq. (38). The path integral over h, a, and µ now factors among the N different sites i, and we obtain Eqs. (36) and (37): Φ n [G, F ] ≡ β 2 2 rr 1 0 dτ dτ G rr (τ, τ ) p G r r (τ , τ ) p − F rr (τ, τ )G r r (τ , τ ) + ln r Dh r Da r Dµ r e −S (eff) [h,a,µ] .
(37) Note that the remaining integration over h, a, and µ is for a single site. In return, the action for that site has couplings between different replicas and times.
Note that for the l = 0 terms, the {IA} sum factors into two separate sums, one for each circle. Furthermore, the sums over k 1 and k 2 are then precisely those that gave us E[Z p ]. Thus For l = 1, η {IA} = 0. This is because eachÔ A I is traceless, and the two operators paired between the circles are each left unpaired in their respective traces.
For l = 2, let us count the powers of N . Assume that k 1 , k 2 ∼ O(1) as well, so that we can ignore minus signs as discussed above. Yet we do still need to ensure that every factor ofσ α i occurs in pairs to survive the trace. This restricts the number of sums over spin indices in {IA} to pk 1 + pk 2 + p: each of the multi-indices within each circle can be summed freely, but the two which connect the circles must have every index paired with each other. The counting for other l ∼ O(1) is analogous, and once we include the number of contractions, we have Note that, at least for p > 2, all l = 0 are suppressed by powers of N relative to l = 0. If we were to naively sum this expression over all (l, k 1 , k 2 ), we would be led to believe that Var[Z p ]/E[Z p ] 2 → 0 as N → ∞, regardless of β. Yet we have proven in the main text that this cannot be true. The resolution, as also discussed in the main text, is that Eq. (B10) holds only for l, k 1 , k 2 ∼ O(1), whereas we need to sum over all values at fixed N . Once (l, k 1 , k 2 ) become comparable to N , not only do anticommuting operators begin to matter, but the combinatorics of the chord diagrams changes. This second point is demonstrated explicitly in the main text.