# A Bosonic Model of Quantum Holography

Brian Swingle  Department of Physics, Brandeis University, Waltham,
Massachusetts 02453, USA  Mike Winer  Joint Quantum Institute, Department of
Physics, University of Maryland, College Park, Maryland 20742, USA

###### Abstract

We analyze a model of qubits which we argue has an emergent quantum
gravitational description similar to the fermionic Sachdev-Ye-Kitaev (SYK)
model. The model we consider is known as the quantum $q$-spin model because it
features $q$-local interactions between qubits. It was previously studied as a
model of a quantum spin glass, and while we find that the model is glassy for
$q=2$, $q=3$, and likely $q=4$, we also find evidence for previously
unexpected SYK-like behavior for the quenched free energy down to the lowest
temperatures for $q\geq 5$. This SYK-like physics includes power-law
correlation functions and an extensive low temperature entropy, so we refer to
the model as Spin SYK. The model is generic in that it includes all possible
$q$-body couplings, lacks most symmetries, and has no spatial structure, so
our results can be construed as establishing a certain ubiquity of quantum
holography in systems dominated by many-body interactions. Furthermore, we
discuss a generalized family of models which includes Spin SYK and which
provably exhibit SYK-like physics in the solvable limit of large local Hilbert
space dimension. We also comment on implications of a bosonic system with SYK-
like properties for the study of holography, Hamiltonian complexity, and
related topics.

##  1 Introduction

In this paper we study the low-energy physics of generic $q$-local quantum
systems. In their simplest form, these are Hamiltonians that include all
possible $q$-body interactions between $N$ sites with minimal additional
structure. The models we consider are bosonic in that their Hilbert space has
a tensor product structure built from the sites such that operators from
different sites commute. Variants of such “mean-field” models have a long
history, including as models of quantum spin glasses [1, 2, 3, 4, 5, 6, 7, 8],
as models for investigating Hamiltonian complexity [9, 10, 11, 12], and, when
generalizing to fermions, as simple models of holography [13, 14, 15, 16, 17,
18, 19, 20, 21, 22, 23]. Most of this paper is concerned with what we call the
Spin SYK model, for which the sites are qubits. We argue that Spin SYK is
glassy at low energies when $q=2,3,4$ but at larger $q$ exhibits holographic
dynamics similar to the fermionic Sachdev-Ye-Kitaev model, hence the moniker.
Even for smaller $q$ for which the physics appears to be glassy at the lowest
energies, Spin SYK can still exhibit SYK-like physics over a window of
energies. Observations of SYK-like physics at intermediate energies has also
been reported in the $q=2$ random Heisenberg model [24, 25].

We also discuss a generalization of Spin SYK, the $q$-local Gauged Clusters
model, in which the qubits are replaced by sites with local Hilbert space
dimension equal to a power of $2$ and with interactions drawn from a large but
proper subset of all local operators. The model is designed to be similar to
random SU($M$) quantum magnets, which are solvable at large $M$ [13]. In the
limit of large local Hilbert space dimension, this $q$-local Gauged Clusters
also becomes solvable and is SYK-like for any $q\geq 2$, although it suffers
from a known low-temperature instability [1, 13, 25] towards glassiness when
$q=2$ and the local Hilbert space is finite.

Our results indicate that the holographic low-energy physics of the Sachdev-
Ye-Kitaev (SYK) model is actually more generic then previously known. Compared
to Spin SYK, the presence of low temperature holographic dynamics appears to
be only slightly more robust in SYK, which is not glassy at $q=4$ [26] and for
which only $q=2$ is non-holographic. The SYK model exemplifies many features
of strong coupling dynamics and provides a simple realization of holographic
duality [27], and these physical properties, combined with its analytical
tractability, have led to its intense study. Compared to SYK, the bosonic
models we consider are interesting for several reasons: (1) they might be
easier to implement in quantum simulations, especially in the near term, (2)
one can potentially get SYK-like physics at intermediate energies even when
$q=2$, and (3) they provide models of holography in which all the microscopic
degrees of freedom are bosonic. Spin SYK is also interesting for its potential
applications in Hamiltonian complexity, e.g. [12], and as a launching point
for further model development. Finally, Spin SYK has the interesting property
that SYK dynamics appear at a quenched saddle point, the one which dominates
the dynamics and thermodynamics for a random couplings. The annealed saddle
point, which dominates the disorder-averaged partition function, is still
glassy. This is in contrast to the SYK model, where the quenched and annealed
saddle points are the same.

As we have said, the bulk of the paper studies the Spin SYK model, which was
introduced in [8] (where it was called a quantum $p$-spin glass model) and
further studied in [28, 29, 30, 31], with Hamiltonian given by

$H=\sum_{r_{1}\mu_{1}\cdots r_{q}\mu_{q}}J_{r_{1}\mu_{1}\cdots
r_{q}\mu_{q}}\sigma_{r_{1}\mu_{1}}\cdots\sigma_{r_{q}\mu_{q}}$ (1)

where $r=1,\cdots,N$ refers to the site and $\mu=1,\cdots,3$ refers to the
spin component (Pauli matrix). The $J$ couplings are Gaussian random variables
with zero mean and variance

$\mathbb{E}\left[J_{r\mu\cdots}^{2}\right]=\frac{(q-1)!J^{2}}{N^{q-1}}.$ (2)

We highlight in particular the recent independent study [30] which considered
a variant of (1) with $q=4$ and only two Pauli components
$\sigma_{x},\sigma_{y}$ instead of all three and gave evidence from finite
size numerics that the physics was SYK-like over a wide range of energies. As
we discuss in Section 9, many of our conclusions will apply to the two
component case as well.

The main assertion of this paper is that model (1) displays power-law decay in
its low temperature correlation functions for sufficiently large $q$. In the
rest of this section, we will define the model and the primary observables of
interest. We will also discuss what we mean when we calculate the disorder-
averaged properties of this model and overview our analysis strategy.

In Section 2 we will write down the action for a path integral formulation of
the Spin SYK thermal partition function in terms of a correlation function $G$
and a self-energy $\Sigma$. We work in the large $N$ limit in which the path
integral can be evaluated using saddle-point methods. We will differentiate
the Spin SYK action to get saddle-point equations for $G$ and $\Sigma$.

In Section 3 we will discuss possible solutions of these equations. We will
prove the existence of two saddle-points of the single-replica action: an
annealed solution with $G\to 1$ at low temperatures which dominates the naive
average of the partition function, and a possible quenched solution which
could control the free energy one would actually measure. This solution is
replica-diagonal and gapless, and obeys a power-law decay.

Section 4 examines the large-$q$ limit of the Spin SYK model, and shows that
it obeys a Liouville equation just like the traditional SYK model. This means
that the large $q$ properties are essentially identical to that of the SYK
model.

Section 5 studies numerical solutions to the single-replica saddle-point
equations for general $q$. For $q=2$, we find only the annealed solution. For
$q>2$, we find the annealed solution and another solution which we refer to as
SYK-like. Moreover, even as low as $q=3$, the SYK-like solution to the
Liouville equations is a good qualitative fit for (an) exact solution. We also
find that the SYK model’s conformal solution is a good fit for the SYK-like
solution at low temperature.

Section 6 goes over exact-diagonalization results, which suggest that the
system is glassy at $q=2,3,4$ and SYK-like for $q\geq 5$. We study a sparse
version of Spin SYK which allows us to address larger $q$ efficiently. Our
main strategy is to measure the Edwards-Anderson order parameter. When it is
zero, then the system should be replica diagonal and the power-law solution
dominates. When it is non-zero, the system is presumably glassy. We also look
at low-energy level spacings, finding results consistent with SYK-like physics
and inconsistent with a spin glass for $q\geq 5$.

Section 7 deals with a generalization of Spin SYK in which the qubits are
replaced by sites of Hilbert space dimension $2^{M/2-1}$ for an even integer
$M\geq 4$. These sites have an interpretation as $M$ Majorana fermions, with
the local fermion parity gauged to produce a bosonic local Hilbert space,
hence $q$-local Gauged Clusters. In this picture, we make a set of all
possible fermion bilinears (which are bosonic operators) and consider generic
$q$-local interactions built from these operators. Spin SYK is the special
case $M=4$. At large $M$, the model is solvable and exhibits SYK-like physics
over a wide temperature range for any $q\geq 2$. We conjecture that this
behavior persists all the way down to $M=4$ for sufficiently large $q$.

In Section 8 we elaborate on just how similar the low energy physics is
between SYK and Spin SYK. We discuss the ground state entropy, the emergence
of spontaneously and explicitly broken time-reparameterizations, the
Schwarzian action, and the possible effects of weakly irrelevant operators.

Finally, in Section 9 we conclude with a brief outlook that highlights many
questions remaining to be addressed. A few appendices contain the derivation
of the Spin SYK action and details of the numerical calculations.

This paper represents the first part of an ongoing investigation. We plan to
report on a detailed study of multi-replica saddles in a future work.

###  1.1 Model and Observables

We consider a model defined on $N$ spin-1/2 degrees of freedom. We also fix an
integer $q>1$, the “locality”, which determines the degree of the interaction.
The terms in the Hamiltonian are obtained by considering all possible $q$-body
subsets of the spins and an assignment of $\sigma_{\mu}$ ($\mu=1,\cdots,3$)
for each spin. This yields an ensemble of Hamiltonians where each instance
takes the form

$H=\sum_{r_{1}\mu_{1}\cdots r_{q}\mu_{q}}J_{r_{1}\mu_{1}\cdots
r_{q}\mu_{q}}\sigma_{r_{1}\mu_{1}}\cdots\sigma_{r_{q}\mu_{q}},$ (3)

with $r=1,\cdots,N$ referring to the spin and $\mu=1,2,3$ referring to the
component. The model does not have spin-rotation symmetry (except in a
statistical sense), and we typically use spin and qubit terminology
interchangeably.

Our goal in this paper is to analyze the Hamiltonian ensemble defined in (3).
Our primary interest is to obtain the energy and entropy in the Gibbs state
$e^{-\beta H}$ as a function of $\beta$, averaged over realizations of the $J$
couplings. We focus on the low temperature regime of $N\to\infty$ and $\beta
J$ large.

Crucial to our considerations will be the Green’s function or correlation
function,

$G(\tau_{1},\tau_{2})=\frac{1}{N}\sum_{r\mu}\mathbb{E}\left[\langle\sigma_{r\mu}(\tau_{1})\sigma_{r\mu}(\tau_{2})\rangle_{\beta}\right],$
(4)

where we take both a thermal and an ensemble average. The arguments of the
Pauli matrices indicate imaginary time values, with $O(\tau)=e^{\tau
H}Oe^{-\tau H}$ for $\tau\in[0,\beta)$. As we will see below, the low-
temperature fate of the model is closely tied to the properties of this
correlation function. In our work, $G(\tau_{1},\tau_{2})$ will be translation
invariant, so that $G(\tau_{1},\tau_{2})=G(\tau_{1}-\tau_{2},0)\equiv
G(\tau_{1}-\tau_{2})$.

At very large $\beta$, an important role is played by the asymptotic value
$G(\infty)$. From the definition of $G$, this is equal to the ensemble average
of the “Edwards-Anderson” order parameter in the ground state,

$G(\infty)=\frac{1}{N}\sum_{r\mu}\langle\text{gs}|\sigma_{r\mu}|\text{gs}\rangle^{2}\equiv
q_{\text{EA}}.$ (5)

This notation is conventional (although our normalization of $q_{\textrm{EA}}$
is slightly different); $q_{\text{EA}}$ has nothing directly to do with $q$,
the locality. A non-zero value of $q_{\text{EA}}$ indicates a frozen pattern
of local expectation values in the ground state, which is a hallmark of a
glass. In contrast, when $q_{\text{EA}}$ is zero at $\beta=\infty$, at least
the equilibrium thermodynamics of the model is not glassy. An important
implication of $q_{\text{EA}}>0$ is that if we have multiple “replicas” of the
system with identical $J$ couplings, then the frozen pattern of expectation
values will have some similarity or “overlap” between different replicas. This
overlap is also controlled by $q_{\text{EA}}$ at large $\beta$.

###  1.2 Quenched vs Annealed

We are interested in the disorder-averaged free energy, $F_{q}$ [32]. This is
a “quenched” average, meaning we compute the free energy for a given sample (a
realization of the $J$ couplings), and then average over samples,

$F_{q}\equiv\int dJP(J)F(J)\equiv\mathbb{E}\left[F\right].$ (6)

Since $F$ is related to the partition function $Z=\tr(e^{-\beta H})$ via

$F=-\frac{1}{\beta}\ln Z,$ (7)

this quenched average corresponds to averaging $\ln Z$,

$F_{q}=-\frac{1}{\beta}\mathbb{E}\left[\ln Z\right].$ (8)

As is well known, it can be challenging to compute the average of $\ln Z$. The
standard method is to use a replica trick in which we introduce $n$ copies or
replicas of the system, compute the average of $Z^{n}$, and then take $n\to 0$
using the identity

$\mathbb{E}\left[\ln Z\right]=\lim_{n\to
0}\frac{\mathbb{E}\left[Z^{n}\right]-1}{n}.$ (9)

Although the replica trick is known or conjectured to have issues in some
cases [33, 34], we will assume it is valid and that we can exchange the order
of limits of $N\to\infty$ and $n\to 0$.

Even with this assumption, evaluating the “replica limit” $n\to 0$ is still
typically challenging. Our approach is to express the partition function as a
path integral. Then the disorder average of $Z^{n}$ is straightforward to
carry out and we are left with a “replicated” path integral to evaluate. In
the large $N$ limit, we write this replicated path integral schematically as

$\mathbb{E}\left[Z^{n}\right]=\int\mathcal{D}\phi e^{-NI_{n}[\phi]}$ (10)

for some general integration variables $\phi$ and “replicated action” $I_{n}$.
Taking the large $N$ limit suggests looking among the saddle points of $I_{n}$
for the physical saddle. However, whereas one might naively look for saddles
of smallest action, there is a well known prescription due to Parisi that
states the correct answer in the replica limit is actually to look for saddles
that are maximal [35], at least with respect to “replica order parameters”.
See Eq. (113) for details of the replicated action for Spin SYK, written in
terms of a variable $G$ which will reduce to the physical correlation function
at the physical saddle point and a corresponding self-energy $\Sigma$, which
jointly play the role of $\phi$.

A far simpler quantity to compute is the logarithm of the average of $Z$. This
order of operations defines the “annealed” free energy,

$F_{a}=-\frac{1}{\beta}\ln\mathbb{E}\left[Z\right],$ (11)

and for SYK we have the remarkable fact that

$\lim_{N\to\infty}\frac{F_{q}}{N}=\lim_{N\to\infty}\frac{F_{a}}{N}.$ (12)

This is an enormous technical simplification in that we only ever have to
consider a single replica, at least for thermodynamic quantities. The way this
works for SYK is that the relevant saddle point of the replicated action,
$\phi_{*}$, has no correlation amongst the replicas. In this case,
$I_{n}(\phi_{*})=nI_{1}(\phi_{*})$ and the replica limit can be taken without
trouble. Thus, the quenched free energy is determined from $I_{1}(\phi_{*})$.
In this situation, we say the relevant saddle is “replica diagonal” (no
correlations between replicas) thus “replica symmetric” (the $S_{n}$ symmetry
among the replicas is unbroken).

Unfortunately, it was recently shown that the quenched and annealed free
energies must differ for a large class of bosonic models below a model-
dependent (but order one) temperature [29]. The Spin SYK model satisfies the
conditions of this theorem, so $F_{q}\neq F_{a}$ at sufficiently low
temperature. In general, for systems like the Spin SYK model, at low
temperatures the annealed free energy will linear in $\beta$, while the
quenched free energy will be the ground state energy, which is independent of
$\beta$. This implies that different saddle points compute the annealed and
quenched free energies. The annealed saddle is given by a function
$G_{a}(\tau)$ which goes to some non-zero constant for large $\tau$, a
constant which is not directly related to $q_{\text{EA}}$. The quenched saddle
$G_{q}(\tau)$ must instead approach $q_{\textrm{EA}}$, the overlap between
replicas. If this $q_{\textrm{EA}}$ is zero, then the quenched solution is
replica diagonal and goes to zero at long times.

This multitude of solutions does not imply that the model is necessarily
glassy, although one often thinks of glassiness when quenched and annealed
energies differ. But it does mean that we must deal with the replicated action
and try to determine the correct saddle point, even if the correct solution is
ultimately replica diagonal. We employ a recently proposed “minimax” procedure
which extends Parisi’s proposal [36]. This procedure calls for first
maximizing $I_{n}$ over replica order parameters (in this case the
correlations between replicas) and then minimizing it over conventional order
parameters (in this case the single-replica correlation).

###  1.3 Analysis Strategy

In principle, the procedure is clear: write down the replicated action, find
its saddle points, and carry out the minimax procedure to determine which
saddle controls the quenched free energy. This calculation is possible but
somewhat involved, so we are deferring it to future work. Here, we follow a
different strategy in which we use finite size numerical results to tell us
whether the quenched saddle is replica diagonal or not. When it is, we can
obtain the quenched free energy from single-replica saddles, which we analyze
in detail.

The key point is that, for disordered systems which are not in a glassy phase,
for example, the Transverse Field Sherrington Kirkpatrick (TFSK) [37, 38]
model at large transverse field, it can be that the minimax saddle is replica
diagonal but is distinct from the annealed saddle point (which is also replica
diagonal). In this case, the minimax saddle can be obtained from the single
replica equations of motion. It will always have higher action compared to the
annealed saddle at sufficiently low temperature, nevertheless, the minimax
procedure instructs us to take the minimax saddle over the annealed saddle
when computing the quenched free energy.

Our strategy thus has two parts. First, we construct an inventory of possible
single-replica saddles. This turns out to be a small set, containing just the
annealed saddle and another SYK-like power-law saddle provided $q>2$. This is
the business of Sections 2,3 4, and 5. Second, we use exact diagonalization to
directly study the overlap $q_{\textrm{EA}}$ in the ground state of sparse
finite $N$ instances of the model. This provides an indirect way to determine
whether the saddle which controls the quenched free energy is replica diagonal
or not, although it relies on an extrapolation to large $N$. This is the
business of Section 6.

If $q_{\textrm{EA}}$ is zero in the thermodynamic limit, then whatever saddle
is chosen by minimax at large $\beta$ must be replica diagonal. Since there is
only one replica diagonal option (remembering that the annealed solution
cannot compute the quenched free energy), we know that it must be the SYK-like
saddle which controls the quenched free-energy. If $q_{\textrm{EA}}$ is non-
zero in the thermodynamic limit, then the model is glassy. Our strategy is
summarized in Figure 1.

Figure 1: The strategy we follow in this paper is in grey. The top flowchart
of blue boxes is the general procedure proposed in [36] in which we construct
the replicated action, carry out the minimax procedure, and then use the
quenched saddle to determine the free energy. While this procedure should be
possible to carry out, here we adopt an alternate procedure in which we try to
determine $q_{\textrm{EA}}$ using an extrapolation from finite size numerics
(orange box). When $q_{\textrm{EA}}=0$, we can restrict to single-replica
saddles in which case our results indicate that the only possibility is the
SYK-like saddle.

The main failure modes of this approach are two-fold. First, we might miss a
relevant single-replica saddle. As we discuss below, we could find no other
solutions and rule out large classes of solutions including exponentially
decaying ones, but we cannot at this time rigorously prove that there are no
additional solutions. Second, the sparse finite $N$ results might be
misleading, either because of the finite sparsity or the finite system size.
For $q\geq 5$, we observed an exponential decrease of $q_{\textrm{EA}}$ with
system size and saw no evidence of significant deviations from this behavior,
but cannot rule out a surprise reversal at large $N$.

##  2 The Action of Spin SYK

We begin our analysis by using spin coherent states [39, 40] to develop a path
integral expression for disorder-average of the partition function
$\mathbb{E}\left[Z\right]$. For the spin-1/2 Pauli matrices in the Spin SYK
Hamiltonian, we have the map $\sigma_{\mu}\rightarrow 3\Omega_{\mu}$ where
$\Omega_{\mu}$ is a unit vector labeling a spin coherent state. We include all
the overlaps of the spin coherent states into a measure
$\mathcal{D}^{N}\Omega$, c.f. (107). The single replica average partition
function is then

$\mathbb{E}\left[Z\right]=\int\mathcal{D}^{N}\Omega\exp\left(\frac{\mathbb{E}\left[J_{r\mu\cdots}^{2}\right]}{2q!}\sum_{A}\left[\int_{0}^{\beta}d\tau
s_{r_{1}\mu_{1}}(\tau)\cdots s_{r_{q}\mu_{q}}(\tau)d\tau\right]^{2}\right),$
(13)

where $s_{rx}(\tau)=3\sin\theta_{r}(\tau)\cos\phi_{r}(\tau)$,
$s_{ry}(\tau)=3\sin\theta_{r}(\tau)\sin\phi_{r}(\tau)$ and
$s_{rz}(\tau)=3\cos\theta_{r}(\tau)$. This sort of path integral was written
down for several closely related models in [29].

We next introduce $G,\Sigma$ variables with the $\Sigma$ variable enforcing
the constraint that

$G(\tau_{1},\tau_{2})=\frac{1}{N}\sum_{r\mu}s_{r\mu}(\tau_{1})s_{r\mu}(\tau_{2}).$
(14)

This will correspond to the correlation function
$\langle\sigma_{x}(\tau_{1})\sigma_{x}(\tau_{2})\rangle+\langle\sigma_{y}(\tau_{1})\sigma_{y}(\tau_{2})\rangle+\langle\sigma_{z}(\tau_{1})\sigma_{z}(\tau_{2})\rangle$,
averaged over sites. Up to $1/N$ corrections, the path integral is

$\mathbb{E}\left[Z\right]=\int\mathcal{D}^{N}\Omega\mathcal{D}\Sigma\mathcal{D}G\exp\left(N\frac{J^{2}}{2q}\int
G^{q}d\tau^{2}-\frac{N}{2}\int\Sigma\left[G-\frac{1}{N}\sum_{r\mu}s_{r\mu}(\tau_{1})s_{r\mu}(\tau_{2})\right]d\tau^{2}\right).$
(15)

The resulting integral is still hard to solve because of the spin bilinear. So
we take one more step of introducing fluctuating magnetic fields
$h_{x,y,z}(\tau)$, one for each spin, each with covariance
$\Sigma(\tau_{1},\tau_{2})$. Then at last the path integral is

$\begin{split}\mathbb{E}\left[Z\right]=\int
D^{N}\Omega\mathcal{D}\Sigma\mathcal{D}G\mathcal{D}^{N}h\exp\left(N\frac{J^{2}}{2q}\int
G^{q}d\tau^{2}-\frac{N}{2}\int\Sigma
Gd\tau^{2}+\int\sum_{r\mu}s_{r\mu}h_{r\mu}\right)=\\\
\int\mathcal{D}\Sigma\mathcal{D}G\exp\left(N\frac{J^{2}}{2q}\int
G^{q}d\tau^{2}-\frac{N}{2}\int\Sigma
Gd\tau^{2}+N\log\langle\mathcal{P}e^{\int\sum_{\mu}\sigma_{\mu}h_{\mu}(\tau)d\tau}\rangle_{\Sigma}\right)\end{split}$
(16)

where $\mathcal{D}h$ is a Gaussian measure with covariance $\Sigma$ and
$\langle\mathcal{P}e^{\int\sum_{\mu}\sigma_{\mu}h_{\mu}(\tau)d\tau}\rangle_{\Sigma}$
denotes a single spin path integral averaged over magnetic fields. For a full
writing of the path integral, see appendix A.

The equations of motion (EOM) or saddle-point equations are

$\begin{split}\Sigma=J^{2}G^{q-1}\,\,\,\,(\text{$G$ EOM})\\\
G=\langle[\sigma_{x}(\tau_{1})\sigma_{x}(\tau_{2})+\sigma_{y}(\tau_{1})\sigma_{y}(\tau_{2})+\sigma_{z}(\tau_{1})\sigma_{z}(\tau_{2})]\rangle_{\Sigma}\,\,\,\,(\text{$\Sigma$
EOM})\end{split}$ (17)

where $\langle\cdots\rangle_{\Sigma}$ denotes a single-spin correlator in the
presence of a random (imaginary) time-dependent field $h(\tau)$ with
covariance $\Sigma$. This correlator should be evaluated as an annealed
average over $h$, meaning we separately average the numerator and denominator.
All the saddles we consider will have time-translation symmetry, meaning
$G(\tau_{1},\tau_{2})=G(\tau_{1}-\tau_{2})$ and
$\Sigma(\tau_{1},\tau_{2})=\Sigma(\tau_{1}-\tau_{2})$, unless otherwise noted.

###  2.1 Sampling Approach

Practically speaking, it is still hard to deal analytically with the
fluctuating $h$ field. However, we can adopt a sampling approach. We divide
$\tau\in[0,\beta)$ up into $N_{\tau}$ intervals of length
$\Delta\tau=\beta/N_{\tau}$. One sample corresponds to $3N_{\tau}$ numbers
$h_{\mu}(j\Delta\tau)$ drawn from a Gaussian distribution with covariance
$\delta_{\mu,\mu^{\prime}}\Sigma(j\Delta\tau,j^{\prime}\Delta\tau)$. We then
construct the single spin path integral

$Z_{1}=\tr\prod_{j}e^{\Delta\tau\sigma_{\mu}h_{\mu}(j\Delta\tau)}.$ (18)

By averaging over $h$, we get $\mathbb{E}\left[Z_{1}\right]_{h}$. We also need
the analog of $Z_{1}$ but with the operators inserted, $Z_{1,\text{op}}$, and
its average over $h$, $\mathbb{E}\left[Z_{1,\text{op}}\right]_{h}$. The
correlator $G$ is

$G=\frac{\mathbb{E}\left[Z_{1,\text{op}}\right]_{h}}{\mathbb{E}\left[Z_{1}\right]_{h}}=\langle\sigma_{\mu}(\tau_{1})\sigma_{\mu}(\tau_{2})\rangle_{\Sigma}.$
(19)

We can approximate both numerator and denominator via sampling. This is an
annealed quantity where we separately average numerator and denominator.

The basic algorithm to solve the EOM starts with a guess for $G$ and $\Sigma$,
uses $\Sigma$ and sampling to estimate the corresponding $G$, updates $G$ and
$\Sigma$ using the other equation of motion and some update rule (e.g. a
linear combination of old and new), and then repeats until convergence. There
are further details in the numerical implementation. For example, the update
rule can be changed dynamically to aid convergence. We will present data
obtained from scheme below in Sec. 5.

##  3 Solving the Saddle-Point Equations

###  3.1 Annealed Solution

We now discuss solutions to the EOM. Before considering more complicated
saddles, let us first discuss the saddle which computes the annealed free
energy, and dominates the disorder-averaged partition function at large
$\beta$. In order to make the annealed free energy proportional to $\beta$ at
large $\beta$, this saddle must have the property that $G(\tau)$ and
$\Sigma(\tau)$ have a non-vanishing limit as $\tau,\beta\to\infty$. If we make
the so-called static approximation, neglecting the time-dependence of $G$ and
$\Sigma$, then we see that the $G$ EOM is solved by

$G=1,\Sigma=J^{2}.$ (20)

At very large $\beta$, these values give the large $\tau$ limit of $G(\tau)$
and $\Sigma(\tau)$.

To see that the $G=1$ solution also satisfies the $\Sigma$ EOM at large
$\beta$, consider the single spin problem with fluctuating field $h$. Since
$\Sigma$ is independent of time (which is approximately valid for most of
$\tau$ at large $\beta$), the fluctuating field $h$ can be taken to be uniform
in imaginary time. A Gaussian distribution with covariance
$\Sigma(\tau,\tau^{\prime})=\text{constant}$ can only generate constant-in-
time $h$ fields.

Thus, we must analyze a single spin in a uniform random field where each
component has variance $J^{2}$ in the limit of large $\beta$. This is
straightforward, as $\beta\to\infty$ simply projects onto the ground state of
the $h$ field. In this ground state, $G$ will be unity since, after a rotation
of the axes to align the $x$-axis with the random field, only $\sigma_{x}$ has
a non-vanishing contribution,

$G\approx\frac{e^{(\beta-\tau)|h|}\langle\sigma_{x}\rangle
e^{\tau|h|}\langle\sigma_{x}\rangle}{e^{\beta|h|}}=\langle\sigma_{x}\rangle^{2}=1.$
(21)

The actual annealed free energy can be obtained via a simple trick. Being
dimensionless, $\mathbb{E}\left[Z\right]$ depends only on the dimensionless
combination $\beta J$. Hence, the thermodynamic identity

$E_{a}=-\partial_{\beta}\ln\mathbb{E}\left[Z\right]$ (22)

can be converted into

$E_{a}=-J\partial_{\beta
J}\ln\mathbb{E}\left[Z\right]=-\frac{J}{\beta}\partial_{J}\ln\mathbb{E}\left[Z\right].$
(23)

Now, dependence on $J$ comes from two sources, the explicit dependence in
front of $G^{q}$ and the implicit dependence via the saddle point solution.
However, the derivative of these fields gives no contribution since the EOM is
obeyed,

$\partial_{J}\text{action}=\text{explicit
term}+\text{EOM}_{G}\partial_{J}G+\text{EOM}_{\Sigma}\partial_{J}\Sigma.$ (24)

Hence, we need only consider the explicit dependence. The result is

$E_{a}=-\frac{J^{2}}{q}\int d\tau G_{a}(\tau)^{q},$ (25)

where again $G_{a}(\tau)\to 1$ as $\tau\to\infty$. Note that this equation is
true beyond the static limit, and analogous identities relate other saddles
with the corresponding energies.

Since $G_{a}$ is approximately constant for most of $\tau$, it follows that

$E_{a}=-\frac{J^{2}}{q}\left(\beta+\cdots\right)$ (26)

at large $\beta$. Because the annealed energy diverges to negative infinity as
$\beta\to\infty$, it cannot possibly be equal to the quenched energy as the
latter is necessarily bounded given the bounded spectrum of the Spin SYK
Hamiltonian. Hence, the annealed free energy, while correctly computing
$\mathbb{E}\left[Z\right]$, fails to capture the quenched physics at
sufficiently large $\beta$.

###  3.2 The Asymptotic Value of $G$

For the annealed solution at large $\beta$, we found $G(\infty)=1$. What other
possible values can $G(\infty)$? Here we argue that the only other possible
value for a solution of the single-replica EOM is $G(\infty)=0$. Note that
other values are possible for solutions of the multi-replica EOM. Even if
other asymptotic values $G(\infty)$ were possible for single-replica
solutions, none of these could control the quenched free energy since they
would all have diverging energy according to (25).

We make use of another way of expressing the energy,

$E=\frac{NG^{\prime}(0^{+})}{4q},$ (27)

which follows from $\sum_{r\mu}[H,\sigma_{r\mu}]\sigma_{r\mu}=4qH$. Both this
formula and (25) compute the energy of a given saddle and agree even when the
saddle is not the dominant one.

We work at large $\beta$ and suppose that $G(\infty)$ is some non-vanishing
constant in this limit. From (25) we learn that $E\sim-J^{2}\beta
G(\infty)^{q}$. From (27) we learn that $G(\tau)$ must vary rapidly near
$\tau=0$ in order to generate such a large $E$, $G^{\prime}(0)\sim-J^{2}\beta
G(\infty)^{q}$. At large $\beta$ and fixed $G(\infty)$, we therefore expect
$G$ to be close to its late time value after a very short time interval of
size $1/(\beta J^{2}G(\infty)^{q})$.

This rapid approach to the late time constant value of $G$ implies a similar
rapid approach of $\Sigma$ to its late time value. Since
$\Sigma-\Sigma(\infty)$ is rapidly decaying, we again appeal to the static
approximation with the expectation that it will be quantitatively accurate at
large $\beta$. The only possible solution with non-vanishing $G(\infty)$ thus
has $G(\infty)=1$. This is the annealed solution we already discussed. The
only other possibility is $G(\infty)=0$, for which the static approximation
must fail.

###  3.3 Non-Annealed Solutions Can’t Decay Exponentially

In this subsection, we argue that the equations of motion don’t have any
rapidly decaying solutions apart from the annealed solution. For the next two
subsections, 3.3 and 3.4, we assume that $G(\infty)=0$ at large $\beta$. The
impossibility of a rapidly decaying which could control the quenched free
energy means that Spin SYK cannot have a gapped “paramagnetic” phase.

As a warmup, we will first show that the model cannot have white noise for
which $\Sigma(\tau_{1},\tau_{2})\propto\delta(\tau_{1}-\tau_{2}$). We start
with the $\Sigma$ equation, written as

$G(\tau_{1},\tau_{2})=\frac{\mathbb{E}\left[Z_{1,\text{op}}\right]_{h}}{\mathbb{E}\left[Z_{1}\right]_{h}}.$
(28)

The denominator can be written as

$\mathbb{E}\left[Z_{1}\right]_{h}=\tr\int\mathcal{D}h_{\mu}\mathcal{P}\exp\left(\int_{0}^{\beta}h_{\mu}(t)\sigma_{\mu}d\tau\right),$
(29)

while the numerator is

$\begin{split}\mathbb{E}\left[Z_{1,\text{op}}\right]_{h}=3\tr\int\mathcal{D}h_{\mu}(0<\tau\leq\tau_{1})\mathcal{D}h_{\mu}(\tau_{1}<\tau\leq\tau_{2})\mathcal{D}h_{\mu}(\tau_{2}<\tau\leq\beta)\\\
\mathcal{P}\exp\left(\int_{0}^{\tau_{1}}h_{\mu}(\tau)\sigma_{\mu}d\tau\right)\exp\left(\int_{\tau_{1}}^{\tau_{2}}\tilde{h}_{\mu}(\tau)\sigma_{\mu}d\tau\right)\exp\left(\int_{\tau_{2}}^{\beta}h_{\mu}(\tau)\sigma_{\mu}d\tau\right),\end{split}$
(30)

with $\tilde{h}$ being $h$ rotated 180 degrees around a given axis. If
$\Sigma$ is a delta function, then
$\int\mathcal{D}h=\int\mathcal{D}\tilde{h}$, so after relabeling our variables
the numerator in equation (30) equals three times the denominator in equation
(29), we find $G(\tau_{1},\tau_{2})=3$. This is obviously incompatible with a
delta function $\Sigma$.

If $\Sigma$ and $G$ decay exponentially, then $h$s at times separated by more
than a few decay times are essentially independent. The boundary effects can
change the integral by at most a multiplicative constant, so we have

$\mathbb{E}\left[Z_{1,\text{op}}\right]_{h}=\alpha\mathbb{E}\left[Z_{1}\right]_{h}$
(31)

for some constant $\alpha$. Hence, starting from our assumption that $G$
decays exponentially, we’ve shown that it plateaus as some fixed $\alpha$,
thus reaching a contradiction and showing that $G$ cannot decay exponentially.

###  3.4 A Power-Law Saddle

We have so far found an annealed solution and just ruled out a solution which
decays exponentially to zero. The remaining possibility is a saddle with
power-law decay.

To explain why such a saddle in fact exists, it is useful analyze the path
integral as expressed in equation (16). From this point of view, the integral
over $\Omega_{\mu}(\tau)$ computes the statistical partition function of a
classical XYZ chain with couplings given by $\Sigma(\tau_{1}-\tau_{2})$. Here
the chain direction is imaginary time and the Green function $G$ is obtained
as the statistical correlation in this XYZ model. The fact that
$\Sigma=J^{2}G^{q-1}$ imposes the non-trivial constraint that the XYZ
couplings ($\Sigma$) are related to the spin correlations ($G$). We also note
that the measure $\mathcal{D}\Omega$ contains strong nearest couplings which
produce perfect ordering in the absence of $\Sigma$, but this is not a self-
consistent solution except at $\beta J=0$.

Let us denote the power-law couplings as

$\Sigma(\tau)\sim\frac{\Sigma_{0}}{|\tau|^{1+s}}$ (32)

for some to-be-determined parameter $s$. This implies a corresponding power-
law in $G$ thanks to the $G$ EOM. To determine $s$ we must obtain another
expression for $G$ by solving the $\Sigma$ EOM, i.e. by solving the
statistical spin chain problem. We work in the large $\beta$ limit for
simplicity. Remember that $\beta$ corresponds to the length of the chain, and
the role of inverse temperature in the stat mech picture of the chain is
played by the magnitude of $\Sigma_{0}$.

As a function of $\Sigma_{0}$, this model has a phase transition from an
ordered state at large $\Sigma_{0}$ to a disordered state at small
$\Sigma_{0}$. The ordered state boasts long-range order in $G$ and so
corresponds to a non-vanishing value of $G$ at large $\tau$. However, this
does not yield a self-consistent solution of the $G$ EOM. At small
$\Sigma_{0}$, the $G$ correlator decays in the same way as $\Sigma$, like
$\tau^{-(1+s)}$ (note the sign of $s$), which is also not a self-consistent
solution. The remaining possibility is that $\Sigma_{0}$ sits right at a phase
transition between the ordered phase and the disordered phase.

Going back to work of Dyson [41], this “long-range fixed point” has been
analyzed extensively (including in more than one dimension) [42, 43, 44, 45,
46, 47]. The most important feature for our purposes is that, at the critical
point, for $\Sigma$ decaying as $\tau^{-(1+s)}$, the corresponding
correlations decay as

$G(\tau)\sim\frac{1}{\tau^{1-s}}.$ (33)

This is why we introduced the $s$ parameter above. Remarkably, the long-range
character of the fixed point means this relation, which is what we would have
for a Gaussian field theory, is valid even with interactions. For $s\leq 0$,
the couplings decay too slowly and the system is always ordered. For $s\geq
1$, the couplings decay too rapidly and one recovers the short-range physics.
There is no phase transition for the short-range XYZ model in 1d. So the
interesting range of $s$ is

$0<s<1.$ (34)

It is also known that for $s<1/2$, the fixed point is Gaussian, while for
$s>1/2$, the Gaussian fixed point is unstable and flows to an interacting
long-range fixed point.

In our case, the $G$ EOM requires

$(1-s)(q-1)=1+s,$ (35)

which yields

$s=1-\frac{2}{q}.$ (36)

The corresponding scalings are

$G\sim\frac{1}{\tau^{2/q}}$ (37)

and

$\Sigma\sim\frac{1}{\tau^{2(1-1/q)}},$ (38)

which are precisely the SYK correlations in the ground state. Of course, this
is essentially a result of dimensional analysis once we know the Gaussian
result for the power-law exponent is robust.

We refer to this saddle as SYK-like. If it is the saddle selected by minimax,
then the physics of Spin SYK at low temperature will share at least some
similarities with the SYK. In particular, it will feature an emergent
conformal symmetry, as the long-range fixed points considered above are known
to enjoy an enhanced conformal symmetry which includes the scaling symmetry
[46].

###  3.5 No $q=2$ SYK-like Saddle

It is interesting to study the SYK-like saddle as a function of $q$. In
particular, for $q=2$ we see that $s=0$. This corresponds to couplings that
decay as $1/\tau$, but such a slow decay actually results in a chain that is
always in an ordered phase in the thermodynamic limit. Hence, we immediately
learn that there actually is no SYK-like saddle for $q=2$. This is good
because the physics of the $q=2$ SYK model is indeed very different. For the
Spin SYK model, it seems the only option for $q=2$ is glassy physics.

For $q=3$, we find $s=1/3$ which is in the Gaussian regime. Hence, we can
replace the XYZ chain by a Gaussian field theory and obtain explicit
Schwinger-Dyson like equations for $G$ in terms of $\Sigma$. For $q=4$, the
value of $s$ is $1/2$, which is right on the boundary between the Gaussian and
interacting regimes. Interestingly, once we have $q=5$ or greater, the fixed
points are all interacting. Hence, there is a significant difference between
the physics of the SYK-like saddle for $q\geq 5$ and $q=3,4$.

###  3.6 Summary of Solutions

We first discussed the annealed solution and then searched for solutions which
could control the quenched free energy. We argued that the only possible
values of $G(\infty)$ were the annealed value and zero. Moreover, any replica
diagonal saddle which controls the quenched free energy must have
$G(\infty)=0$ (otherwise the energy (25) would diverge with $\beta$). Among
replica diagonal solutions with $G(\infty)=0$, we showed that an exponentially
decaying solution was not possible. However, for $q>2$, there is an
interesting solution with power-law correlations.

Hence, we have the following structure for the quenched free energy. For
$q=2$, the only available solution is a replica non-diagonal one (which we did
not construct). This strongly suggests that $q=2$ is generically glassy, and
this is confirmed by evidence from exact diagonalization. This conclusion is
broadly consistent with the presence of glassy physics in many quantum spin
models with random $2$-local couplings, e.g. [13].

For $q>2$, we do have the SYK-like saddle, but we may also have replica non-
diagonal saddles. The correct physics is determined by minimax. We have not
yet carried out the full minimax procedure, so here we rely on evidence from
exact diagonalization (see Sec. 6) to determine the fate of the system for
$q>2$.

To substantiate the general arguments in this section, we first discuss the
single-replica solutions in more detail using an analytical large $q$ approach
(Section 4) and a numerical approach (Section 5). Then we return to the exact
diagonalization analysis in Section 6.

##  4 Large $q$ Limit and the Liouville equation

To gain further insight into the problem (and because it will be useful for
numerically solving the EOM), we can study the large $q$ limit [15, 48]. We
start with the equations of motion

$\Sigma(\tau_{1},\tau_{2})=J^{2}G^{q-1}(\tau_{1},\tau_{2})$ (39)
$G(\tau_{1},\tau_{2})=\langle\sigma_{\mu}(\tau_{1})\sigma_{\mu}(\tau_{2})\rangle_{\Sigma}$
(40)

and we make the ansatz

$\begin{split}G=3\left(1+\frac{g}{q}+\cdots\right),\\\
\Sigma(\tau)=3^{q-1}J^{2}e^{g(\tau)}+\cdots.\end{split}$ (41)

To close the EOM, we need to compute $G$ in terms of $\Sigma$.

This can be done conveniently starting with a formula for the derivatives of
$G$,

$\partial_{\tau_{1}}\partial_{\tau_{2}}G=\frac{\int\mathcal{D}h\tr(\mathcal{P}\left\\{e^{\int
h\cdot\sigma}[h_{\nu}(\tau_{1})\sigma_{\nu}(\tau_{1}),\sigma_{\mu}(\tau_{1})][h_{\rho}(\tau_{2})\sigma_{\rho}(\tau_{2}),\sigma_{\mu}(\tau_{2})]\right\\})}{\int\mathcal{D}h\tr(\mathcal{P}\left\\{e^{\int
h\cdot\sigma}\right\\})}.$ (42)

In the perturbative limit where $G$ is close to $3$, we expand (42)
perturbatively in $\Sigma$. Since the numerator already contains a factor
$\Sigma$ that arises from the $h_{\nu}(\tau_{1})h_{\rho}(\tau_{2})$ factors,
we can neglect the fields in the exponents to obtain

$\partial_{\tau_{1}}\partial_{\tau_{2}}G\approx\frac{\Sigma(\tau_{1},\tau_{2})\tr([\sigma_{\nu},\sigma_{\mu}][\sigma_{\nu},\sigma_{\mu}])}{\tr(1)}=-24\Sigma(\tau_{1},\tau_{2}).$
(43)

On the other hand, we have

$\partial_{\tau_{1}}\partial_{\tau_{2}}G\approx\frac{3}{q}\partial_{\tau_{1}}\partial_{\tau_{2}}g$
(44)

from the large $q$ ansatz, so we can now close the EOM,

$\partial_{\tau_{1}}\partial_{\tau_{2}}g=-8q3^{q-1}J^{2}e^{g}.$ (45)

We introduce the new coupling

$\mathcal{J}^{2}=4q3^{q-1}J^{2}.$ (46)

and make the further ansatz $g(\tau_{1},\tau_{2})=g(\tau_{1}-\tau_{2})$, to
finally obtain

$\partial_{\tau}^{2}g=2\mathcal{J}^{2}e^{g},$ (47)

which is the Liouville equation. The standard solution is

$e^{g}=\left[\frac{\cos\frac{\pi v}{2}}{\cos\left[\pi
v\left(\frac{1}{2}-\frac{\tau}{\beta}\right)\right]}\right]^{2}$ (48)
$\beta\mathcal{J}=\frac{\pi v}{\cos\frac{\pi v}{2}}.$ (49)

Note that this has the right boundary values, with $g$ vanishing at
$\tau=0,\beta$.

Using (27) for the energy, we can immediately obtain the thermodynamics from
the derivative of $g$ at $\tau=0^{+}$. As an example, in the $\beta\to\infty$
limit, the standard solution limits to

$e^{g}=\frac{1}{(\mathcal{J}\tau+1)^{2}}.$ (50)

Near $\tau=0$ we thus have $g=-2\mathcal{J}\tau+\cdots$, and so

$-\frac{E}{N}=-\frac{G^{\prime}(0^{+})}{4q}=\frac{3\mathcal{J}}{2q^{2}}.$ (51)

We also get a large-$q$ expression for the ground state entropy analogous to
the SYK result:

$S_{0}/N\approx\log 2-\frac{3\pi^{2}}{8q^{2}}$ (52)

It is also useful to perform an approximate resummation of the large $q$
expansion in which we replace

$G=3\left(1+\frac{g}{q}+\cdots\right)\to 3e^{g/q}.$ (53)

For example, at large $\beta$ this gives

$G\to\frac{3}{(\mathcal{J}\tau+1)^{2/q}},$ (54)

which clearly displays a scaling behavior with $G\sim\tau^{-2/q}$. This
resummation can be understood as arising from an estimate of (42),

$\partial_{\tau_{1}}\partial_{\tau_{2}}G\approx-8\Sigma(\tau_{1},\tau_{2})G(\tau_{1},\tau_{2}).$
(55)

Using the short-time form of $G$ gives the previous result of $-24\Sigma$, but
more generally we can make the ansatz $G=3e^{g/q}$ and compute

$\partial_{\tau_{1}}\partial_{\tau_{2}}G=\left(\frac{1}{q}\partial_{\tau_{1}}\partial_{\tau_{2}}g+\frac{1}{q^{2}}\partial_{\tau_{1}}g\partial_{\tau_{2}}g\right)G.$
(56)

At large $q$, we neglect the second term and again obtain the Liouville
equation.

##  5 Numerical Solution of the Equations of Motion

Here we discuss approximate numerical solutions of the EOM obtained from an
iterative approach. A similar approach was used many years ago in the $q=2$
case [49]. The precise numerical details are described in App. B. We
considered two kinds of initializations for the iterative solver. The first is
the $J=0$ solution, $G(\tau)=3$, which often converges to the annealed saddle
at sufficiently large $\beta$. For the second, we found it very useful to
initialize with the analytic large-$q$ form obtained above. Specifically, we
exponentiate the standard solution to give

$G_{\text{large}}=3\left[\frac{\cos\frac{\pi v}{2}}{\cos\left[\pi
v\left(\frac{1}{2}-\frac{\tau}{\beta}\right)\right]}\right]^{2/q}.$ (57)

This second type of initialization typically converges to something consistent
with the interesting quenched saddle discussed above, even if the initial
guess is not particularly close, i.e. $q$ is not large.

The basic structure of the solver is as follows. The imaginary time circle is
discretized into $N_{\tau}$ units of length $\Delta\tau=\beta/N_{\tau}$ and we
assign a value of $G$ for each time point $\tau_{\ell}=\ell\Delta\tau$.
Similarly, $\Sigma$ is an $N_{\tau}\times N_{\tau}$ matrix. We enforce
translation symmetry in $G$ and $\Sigma$ and the
$\tau\leftrightarrow\beta-\tau$ symmetry in $G$. The method works by proposing
an update, $G_{\text{prop}},\Sigma_{\text{prop}}$, based on the current
values, $G_{\text{curr}},\Sigma_{\text{curr}}$. The proposed new $\Sigma$ is
simply $\Sigma_{\text{prop}}=J^{2}G_{\text{curr}}^{q-1}$. The proposed new $G$
is obtained by sampling the single-spin problem with magnetic fields of
covariance $\Sigma_{\text{curr}}$. Then $G$ is updated to some linear
combination of $G_{\text{curr}}$ and $G_{\text{prop}}$ and similarly for
$\Sigma$. To assess the degree of convergence, we monitor an error given by

$\varepsilon=\sqrt{\frac{1}{N_{\tau}}\sum_{\ell}(G_{\text{curr}}(\tau_{\ell})-G_{\text{prop}}(\tau_{\ell}))^{2}}+\sqrt{\frac{1}{N_{\tau}}\sum_{\ell}(G_{\text{curr}}(\tau_{\ell})^{q-1}-G_{\text{prop}}(\tau_{\ell})^{q-1})^{2}}.$
(58)

Note that this error is intensive in $N_{\tau}$, so it represents an error per
site. When it is small, then the proposed update is close to the current value
on a site-by-site basis. Practically speaking, we can never make this error
smaller than our sampling error.

In addition to the large $q$ analytics, we also compare with the SYK conformal
solution,

$G_{\text{conf}}=b\left(\frac{\pi}{\beta\sin\frac{\pi\tau}{\beta}}\right)^{2/q},$
(59)

with
$J^{2}b^{q}=\frac{\left(\frac{1}{2}-\frac{1}{q}\right)\tan\frac{\pi}{q}}{\pi}$.
A priori, even if the Spin SYK model has a power-law solution, it would not
have to be exactly equal to the SYK solution, but in practice we find that
this is the case.

Figure 2 and Figure 3 show annealed-type solution for $q=2$ and $q=3$,
respectively. These solutions can be obtained with modest computational
effort. We also see that the static approximation looks good, and we expect it
will only improve with increasing $\beta$. For $q=2$, this is the only type of
solution we have been able to find.

Figure 2: Solution to the EOM with $q=2$, $\beta=4$, $J=1$, $N_{\tau}=80$, and
$N_{s}=1000$ samples. The red dashed line is the initialization (large $q$),
the black line shows the final converged result, and the blue shaded region
gives $3/\sqrt{N_{s}}$ around the black line. The result is consistent with
the annealed saddle in which $G$ approaches unity at long time.  Figure 3:
Solution to the EOM with $q=3$, $\beta=4$, $J=1$, $N_{\tau}=80$, and
$N_{s}=1000$ samples. The red dashed line is the initialization ($J=0$), the
black line shows the final converged result, and the blue shaded region gives
$3/\sqrt{N_{s}}$ around the black line. The result is consistent with the
annealed saddle in which $G$ approaches unity at long time.

For $q=3$, one can find other solution types using a large $q$ initialization.
Such a solution is shown in Figure 4 where we also compare it with the SYK
conformal result. The conformal result is an excellent approximation except
very close to $\tau=0,\beta$ (where the large $q$ ansatz is closer). We also
find that the computational resources need to be significantly increased to
find such solutions reliably, including more time points and more samples to
estimate $G$ from $\Sigma$. Readers familiar with SYK may be surprised that
$\beta=4$ seems already solidly in the conformal regime, but this can be
understood because the energy scales of Spin SYK are enhanced by a large
factor arising from the different normalizations of the operators and the
couplings.

Figure 4: Solution to the EOM with $q=3$, $\beta=4$, $J=1$, $N_{\tau}=300$,
and $N_{s}=10000$ samples. Compared to the previous figures, we have now
dramatically increased the number of time points and the number of samples.
The red dashed line is the initialization (large $q$), the black line shows
the final converged result, and the blue shaded region gives $3/\sqrt{N_{s}}$
around the black line. The blue dashed line is the SYK conformal result
(capped at $3$). We see quite good agreement with the large $q$ result at
short times and the SYK result at long times. Note how the solution goes well
below the annealed value at large $\tau$.

We can find similar solutions for $q=4$ as shown in Figure 5. Here we set
$\beta=1$ instead of $\beta=4$ in the $q=3$ and still find excellent agreement
with the conformal ansatz. This is again demonstrating the very rapid rise in
energy scale with increasing $q$, c.f. (46). We can presumably find such a
solution for any $q$, although the computational cost seems to continue to
increase with increasing $q$, at least with our current approach.

Figure 5: Solution to the EOM with $q=4$, $\beta=1$, $J=1$, $N_{\tau}=200$,
and $N_{s}=10000$ samples. The red dashed line is the initialization (large
$q$), the black line shows the final converged result, and the blue shaded
region gives $3/\sqrt{N_{s}}$ around the black line. The blue dashed line is
the SYK conformal result (capped at $3$). This result is roughly consistent
with the large $q$ form for all $\tau$. We again see good agreement with the
large $q$ result at short times and the SYK result at long times.

Finally, we show an example with $q=5$ and even smaller $\beta=.1$ in Figure
6. For this choice of parameters, we are not in the conformal regime but the
large $q$ ansatz gives a good account of the solution.

Figure 6: Solution to the EOM with $q=5$, $\beta=.1$, $J=1$, $N_{\tau}=200$,
and $N_{s}=10000$ samples. The red dashed line is the initialization (large
$q$), the black line shows the final converged result, and the blue shaded
region gives $3/\sqrt{N_{s}}$ around the black line. The blue dashed line is
the SYK conformal result (capped at $3$). The blue dashed line is the SYK
conformal result (capped at $3$). We have significantly reduced $\beta$
relative to the previous figures to demonstrate the good agreement with the
large $q$ result.

##  6 Finite Size Numerics on Sparse Spin SYK

In this section we define and study a sparse version of Spin SYK using sparse
matrix techniques [23, 50, 51]. We consider a simple variant of the sparse
model defined by randomly pruning all but $\kappa N$ of the ${N\choose
q}3^{q}$ terms from the Hamiltonian. The remaining couplings are Gaussian with
a renormalized variance adjusted so that $\mathbb{E}\left[\tr(H^{2})\right]$
has the same value in the dense and sparse models. For the dense model, the
second moment of $H$ is

$\mathbb{E}\left[\tr(H^{2})/2^{n}\right]=\mathbb{E}\left[J_{r\mu\cdots}^{2}\right]\binom{N}{q}3^{q}=NJ^{2}\frac{3^{q}}{q}.$
(60)

The couplings in the sparse model have variance $\tilde{J}^{2}$ obeying
$\kappa N\tilde{J}^{2}=NJ^{2}\frac{3^{q}}{q}$, which gives

$\tilde{J}^{2}=J^{2}\frac{3^{q}}{\kappa q}.$ (61)

We have observed little variation in observables as $\kappa$ is varied. In the
first set of numerical data discussed below, we set $\kappa=N$, which
corresponds to a model that is still relatively sparse but which at large $N$
is likely dense enough largely replicate the physics of the dense model, e.g.
$1/\kappa$ corrections to intensive quantities should vanish at large $N$.

We consider three observables in the sparse model. 1) The ground state
Edwards-Anderson order parameter [52, 53, 54]:
$q_{\textrm{EA}}=\frac{1}{N}\sum_{r\mu}\left(\bra{\psi}\sigma_{r\mu}\ket{\psi}\right)^{2}$
for the many-body ground state $\psi$. This order parameter serves as a
diagnostic for glassy behavior, being nonzero for glasses but zero for spin
liquids or gapped paramagnets. 2) The gap between the ground state and the
first excited state. This would have three different behaviors for the three
possible phases. For a gapped paramagnet, we would most likely see a sharp
peak. For a spin glass, we would see a continuous distribution down to zero
[55, 56]. For an SYK-like system with both level repulsion [57, 58, 59, 60,
61, 62, 63, 64] and an extensive ground state entropy, one would expect a
smooth distribution going to zero as the gap goes to zero. 3) The ground state
energy. This is a fairly self-explanatory quantity. It provides a sanity check
on our overall understanding of the system and can be used to compare with
large $q$ estimates.

We discuss $q=3$, $q=4$, and $q=5$ in turn, starting with $q=3$. The EA order
parameter is shown in Figure 7 where we see an exponential rise as a function
of system size. This indicates that the $q=3$ ground state is a glass. This
conclusion is supported by the gap histogram in Figure 8. This is exactly what
we would expect if the ground state and the first excited state are distinct
TAP states with a macroscopic barrier to tunneling between them.

Figure 7: The Edwards-Anderson order parameter in the ground state (log scale)
for $q=3$ as a function of system size. Each point is averaged over $2000$
samples. The rapid growth with system size strongly indicates glassiness. The
red line is a linear fit to $\log q_{\textrm{EA}}$. Figure 8: A histogram of
$2000$ samples of the finite size energy gap between the ground state and
first excited state for $q=3$ and $n=14$. The exponential distribution of the
gap and the corresponding lack of level repulsion are consistent with the two
lowest lying states being quantum TAP states.

Next, let us consider $q=4$. This case is more complicated because the model
has an anti-unitary time-reversal symmetry which takes
$\sigma_{r\mu}\to-\sigma_{r\mu}$. For odd $N$, the states organize into
degenerate doublets related by time-reversal, while for even $N$ there is no
exact degeneracy. We generalize the $q_{\textrm{EA}}$ order parameter for even
$q$ to take into account the lowest two states,

$q_{EA}=\frac{1}{2N}\sum_{r\mu}\sum_{ab=1,2}|\langle
a|\sigma_{r\mu}|b\rangle|^{2}.$ (62)

The modified EA order parameter is shown in Figure 9 where we see a slow
decrease as a function of system size. This behavior is suggestive of a non-
glassy ground state, although this case seems marginal compared to the robust
decrease seen below for $q=5$. A non-glassy ground state is supported by the
gap histogram in Figure 10 where we see level repulsion. Nevertheless, the
asymptotic limit at large $N$ is less definite than for $q=3$ and $q=5$. We
also show an extrapolation of the ground state energy in Figure 11. The
proximity to the large $q$ result is synergistic with (27) and our observation
that the large $q$ ansatz works well near $\tau=0$.

Figure 9: The Edwards-Anderson order parameter in the ground state (log scale)
for $q=4$ as a function of system size. Each point is averaged over $500$
samples. We slow decrease with system size which is consistent with the
absence of glassiness. We also see a very clear even-odd effect consistent
with the anti-unitary time-reversal symmetry. The red line is a linear fit to
$\log q_{\textrm{EA}}$. Figure 10: A histogram of $500$ samples of the finite
size energy gap between the ground state and first excited state for $q=4$ and
$n=15$. The clear level repulsion is consistent with a non-glassy low energy
state. Figure 11: The energy per spin plotted versus $1/N$ for $q=4$. The red
line is a linear extrapolation to $1/N=0$ and the blue star is the large $q$
prediction. Each point is $500$ samples.

Finally, let us consider $q=5$. The EA order parameter is shown in Figure 12
where we see an exponential decrease as a function of system size. This
indicates that the $q=5$ ground state is not a glass. This conclusion is
supported by the gap histogram in Figure 13. This is exactly what we would
expect if the ground state and the first excited state had random-matrix-like
level repulsion. We also show an extrapolation of the ground state energy in
Figure 14. The proximity to the large $q$ result is synergistic with (27) and
our observation that the large $q$ ansatz works well near $\tau=0$.

Figure 12: The Edwards-Anderson order parameter in the ground state (log
scale) for $q=5$ as a function of system size. Each point is averaged over
$2000$ samples. The rapid decay with system size strongly indicates the
absence of glassiness. The red line is a linear fit to $\log q_{\textrm{EA}}$.
Figure 13: A histogram of $2000$ samples of the finite size energy gap between
the ground state and first excited state for $q=5$ and $n=14$. The clear level
repulsion is consistent with a non-glassy low energy state. Figure 14: The
energy per spin plotted versus $1/N$ for $q=5$. The red line is a linear
extrapolation to $1/N=0$ and the blue star is the large $q$ prediction. Each
point is $2000$ samples.

To shed further light on the physics, and in particular to try to understand
the fate of $q=4$ at large $N$, we carried out a second computation at
somewhat larger system size, up to $N=17$, with $100$ samples per point. For
these data we set the sparsity parameter to $\kappa=16$. The results for the
EA order parameter for $q=3,4,5,6,7$ are shown in Figures 15, 16, 17, 18, and
19. At these somewhat larger sizes, we see indications of saturation for
$q=3,4$ and robust exponential decays with system size for $q\geq 5$.

Figure 15: The Edwards-Anderson order parameter in the ground state (log
scale) for $q=3$ as a function of system size from $N=13$ to $N=17$. Each
point is averaged over $100$ samples. The approximate saturation with system
size indicates a glassy ground state at large $N$. The red line is a linear
fit to $\log q_{\textrm{EA}}$. Figure 16: The Edwards-Anderson order parameter
in the ground state (log scale) for $q=4$ as a function of system size from
$N=13$ to $N=17$. Each point is averaged over $100$ samples. The approximate
saturation with system size indicates a glassy ground state at large $N$. The
red line is a linear fit to $\log q_{\textrm{EA}}$. Figure 17: The Edwards-
Anderson order parameter in the ground state (log scale) for $q=5$ as a
function of system size from $N=13$ to $N=17$. Each point is averaged over
$100$ samples. The rapid decrease with system size indicates a non-glassy
ground state at large $N$. The red line is a linear fit to $\log
q_{\textrm{EA}}$. Figure 18: The Edwards-Anderson order parameter in the
ground state (log scale) for $q=6$ as a function of system size from $N=13$ to
$N=17$. Each point is averaged over $100$ samples. The rapid decrease with
system size indicates a non-glassy ground state at large $N$. The red line is
a linear fit to $\log q_{\textrm{EA}}$. Figure 19: The Edwards-Anderson order
parameter in the ground state (log scale) for $q=7$ as a function of system
size from $N=13$ to $N=17$. Each point is averaged over $100$ samples. The
rapid decrease with system size indicates a non-glassy ground state at large
$N$. The red line is a linear fit to $\log q_{\textrm{EA}}$.

Our conclusion from this section is that $q=3$ is glassy, $q=4$ is likely
glassy, and $q\geq 5$ shows no sign of glassiness. Although we did not study
$q>5$ in detail, based on all the accumulated evidence and trends, we
conjecture that $q\geq 5$ are all non-glassy and SYK-like down to the lowest
temperatures. For $q=3,4$, the relatively small values of $q_{\textrm{EA}}$
suggest the existence of a wide intermediate energy window where the physics
could still be SYK-like even if the system is glassy at asymptotically low
temperature.

##  7 The $q$-local Gauged Clusters Model

To gain additional analytical insight into the physics, this section studies a
generalization of Spin SYK that we call the $q$-local Gauged Clusters model.
This model is indexed by an even integer $M=4,6,\cdots$ such that $M=4$ is
Spin SYK and the model is solvable at large $M$. The model is very analogous
to the SU($M$) magnets studied by Sachdev and Ye [13], suitably modified since
spin-rotation symmetry plays no role for us. We emphasize that the model is
bosonic for any $M$, but it is convenient to construct it using fermions. We
consider $N$ sites each of which contains $M$ Majorana fermions. These
fermions are denoted $\chi_{r\alpha}$, with $r=1,\cdots,N$ and
$\alpha=1,\cdots,M$, and they obey

$\\{\chi_{r\alpha},\chi_{r^{\prime}\alpha^{\prime}}\\}=\delta_{rr^{\prime}}\delta_{\alpha\alpha^{\prime}}.$
(63)

For each site, we have a fermion parity operator,

$(-1)^{F_{r}}=(-2i)^{M/2}\chi_{r1}\cdots\chi_{rM}.$ (64)

The physical Hilbert space is the even parity sector,

$(-1)^{F_{r}}=1,$ (65)

for each site $r$. The local Hilbert space dimension is thus $d=2^{M/2-1}$.

From each site we can construct fermion bilinears,

$O_{r,\alpha\beta}=-2i\chi_{r\alpha}\chi_{r\beta},$ (66)

which are bosonic operators that commute with all the $(-1)^{F_{r}}$s and
square to the identity. The Hamiltonan is built from these operators as

$H=\sum_{r_{1}\alpha_{1}\beta_{1}\cdots
r_{q}\alpha_{q}\beta_{q}}J_{r_{1}\alpha_{1}\beta_{1}\cdots
r_{q}\alpha_{q}\beta_{q}}O_{r_{1},\alpha_{1}\beta_{1}}\cdots
O_{r_{q},\alpha_{q}\beta_{q}}.$ (67)

The $J$ couplings have variance

$\mathbb{E}\left[J_{\cdots}^{2}\right]=\frac{(q-1)!J^{2}}{N^{q-1}}\left[\frac{\binom{M}{2}^{q}2^{2q+1}}{M}\right]^{-1}.$
(68)

The bizarre numerical factor in square brackets is chosen to make the
equations of motion simple.

In the case $M=4$, one may check that the $O_{r,\alpha\beta}$, of which there
are six per site, can be identified with the Pauli matrices within the two-
dimensional even parity subspace. Hence, apart from the fact that each Pauli
appears twice in the sum over $\alpha<\beta$ (e.g. $O_{12}=O_{34}$ on the
physical Hilbert space), we obtain exactly the Hamiltonian of Spin SYK when
$M=4$. This doubling of the two-fermions operators is a consequence of the
fact that $M=4$, and the four-fermi operator $(-1)^{F_{r}}$ is one on the
physical Hilbert space.

We now study the model at large $M$ and show that it is solvable using
standard SYK technology. In fact, this model is a generalization of coupled
cluster models (e.g. [65]) in which we allow for $q$-body interactions between
the clusters/sites and we gauge the fermion parity on each cluster/site. The
result is a purely bosonic model in the sense discussed above: the total
Hilbert space is a tensor product of local Hilbert spaces and operators on
different tensor factors commute. The fermions $\chi_{\alpha}$ are analogous
to “partons” or “spinons” in the condensed matter language, and they offer a
convenient way to describe the operators of interest, the $O_{\alpha\beta}$.

We discuss the single replica path integral for $q$-local Gauged Clusters in
two stages, first for general $M$ and then at large $M$. Using fermion
coherent states, the path integral for general $M$ is

$\mathbb{E}\left[Z\right]=\int\mathcal{D}^{NM}\chi\left\\{\prod_{r}\frac{1+(-1)^{F_{r}}}{2}\right\\}\exp\left(-\sum_{r\alpha}\frac{1}{2}\int
d\tau\chi_{r\alpha}\partial_{\tau}\chi_{r\alpha}+\frac{J^{2}}{2qN^{q-1}}\frac{M}{\binom{M}{2}^{q}2^{2q+1}}\sum_{r_{1},\alpha_{1}<\beta_{1}\cdots}\left[\int
d\tau O_{r_{1},\alpha_{1}\beta_{1}}\cdots\right]^{2}\right).$ (69)

The factor in curly braces is a projector which restricts to even fermion
parity on each site. We then decouple the interaction term by introducing a
Lagrange multiplier term,

$\frac{MN}{16}\int d\tau
d\tau^{\prime}\Sigma(\tau,\tau^{\prime})\left[G(\tau,\tau^{\prime})-\frac{1}{N\binom{M}{2}}\sum_{r,\alpha<\beta}O_{r,\alpha\beta}(\tau)O_{r,\alpha\beta}(\tau^{\prime})\right],$
(70)

which allows the path integral to be written as

$\mathbb{E}\left[Z\right]=\int\mathcal{D}G\mathcal{D}\Sigma e^{-NI}$ (71)

with

$I=\int d\tau
d\tau^{\prime}\left[-\frac{J^{2}}{2q}\frac{M}{2^{2q+1}}G(\tau,\tau^{\prime})^{q}+\frac{M}{16}\Sigma(\tau,\tau^{\prime})G(\tau,\tau^{\prime})\right]+I_{\text{site}}$
(72)

and

$I_{\text{site}}=-\log\left[\int\mathcal{D}^{M}\chi\left\\{\frac{1+(-1)^{F}}{2}\right\\}\exp\left(-\sum_{\alpha}\frac{1}{2}\int
d\tau\chi_{\alpha}\partial_{\tau}\chi_{\alpha}+\frac{M}{16\binom{M}{2}}\sum_{\alpha<\beta}\int
d\tau
d\tau^{\prime}\Sigma(\tau,\tau^{\prime})O_{\alpha\beta}(\tau)O_{\alpha\beta}(\tau^{\prime})\right)\right].$
(73)

Before proceeding to the large $M$ limit, let’s understand the effect of the
fermion parity projector. The effect of the $(-1)^{F}$ insertion is to change
the boundary conditions of the fermions around the imaginary time circle from
anti-periodic (AP) to periodic (P). Hence, the full single site path integral
is a sum of two terms,

$e^{-I_{\text{site}}}=\frac{Z_{\text{AP}}+Z_{\text{P}}}{2},$ (74)

where $Z_{\text{AP}}$ is a conventional term without $(-1)^{F}$ and
$Z_{\text{AP}}$ is a twisted term $Z_{1,P}$ with $(-1)^{F}$. In the large $M$
limit, this is not expected to change the physics of the model significantly,
meaning that the free energy per fermion should be the same in the even and
odd parity sectors, up to $1/M$ corrections. We first describe the large $M$
limit for the untwisted contribution, which counts both fermion parity
sectors. Then we discuss how the twisted contribution is obtained.

The point of the large $M$ limit is that it allows us to evaluate the path
integral over $\chi_{\alpha}$ in the definition of $I_{\text{site}}$ by saddle
point. We introduce another Lagrange multiplier term,

$\frac{M}{2}\int d\tau
d\tau^{\prime}\Sigma_{\chi}(\tau,\tau^{\prime})\left[G_{\chi}(\tau,\tau^{\prime})-\frac{1}{M}\sum_{\alpha}\chi_{\alpha}(\tau)\chi_{\alpha}(\tau^{\prime})\right],$
(75)

so the single site path integral becomes

$Z_{\text{AP}}=\int\mathcal{D}G_{\chi}\mathcal{D}\Sigma_{\chi}\exp\left(M\log\text{Pf}(\partial_{\tau}-\Sigma_{\chi})+\frac{M}{2}\int
d\tau
d\tau^{\prime}\left[\frac{1}{2}\Sigma(\tau,\tau^{\prime})G_{\chi}(\tau,\tau^{\prime})^{2}-\Sigma_{\chi}(\tau,\tau^{\prime})G_{\chi}(\tau,\tau^{\prime})\right]\right).$
(76)

The EOM for $\Sigma_{\chi}$ and $G_{\chi}$ are, respectively,

$G_{\chi}=\frac{1}{\partial_{\tau}-\Sigma_{\chi}}$ (77)

and

$\Sigma_{\chi}=\Sigma G_{\chi}.$ (78)

The EOM for $G$ and $\Sigma$ are, respectively,

$\Sigma=\frac{J^{2}}{2^{2q-2}}G^{q-1}=J^{2}\left(\frac{G}{4}\right)^{q-1}$
(79)

and

$G=4G_{\chi}^{2}.$ (80)

These equations can be reduced to a closed set involving just $G_{\chi}$ and
$\Sigma_{\chi}$,

$G_{\chi}=(\partial_{\tau}-\Sigma_{\chi})^{-1},$ (81)
$\Sigma_{\chi}=J^{2}G_{\chi}^{2q-1}.$ (82)

From this we learn that $G_{\chi}$ is the correlator of an SYK model with
$q_{\chi}=2q$.

What about the twisted contribution, $Z_{\text{P}}$? Starting with a solution
$G(\tau,\tau^{\prime})$ of the untwisted path integral, we can obtain a
solution of the twisted EOM by modifying $G$ such that it picks up an extra
minus sign whenever $\tau$ or $\tau^{\prime}$ move past $\beta$. The reader
may worry that this choice of time breaks time-translation symmetry, but it is
actually a gauge artifact and time-translation symmetry is preserved for any
gauge invariant observable. This twisted $G$ is
$\tilde{G}(\tau,\tau^{\prime})=Q(\tau)G(\tau,\tau^{\prime})Q(\tau^{\prime})$
for some $Q$ that implements the twist, e.g. $Q=\text{sgn}(\beta-\tau)$ for
$\tau$ near $\beta$. Since this is effectively a pre- and post-multiplication
of $G$ by a diagonal matrix, $Q$, that squares to the identity, it follows
that

$\text{Pf}(\tilde{G}^{-1})=\text{Pf}(QG^{-1}Q)=\text{Pf}(G^{-1})$ (83)

and

$\int\tilde{G}^{q}=\int G^{q},\int\tilde{G}\tilde{\Sigma}=\int G\Sigma.$ (84)

Hence, the value of the saddle point action will be identical. Now, it may
still happen that the full twisted contribution vanishes, but we at least
learn that $Z_{\text{P}}=cZ_{\text{AP}}$ for some value $c$ (possibly zero)
which is only polynomially large in $M$ (and hence cannot change the leading
order free energy per fermion at large $M$).

What all this means is simply that, at large $M$, we can ignore the fermion
parity projector. We can therefore simply stick with the standard anti-
periodic contribution to $e^{-I_{\text{site}}}$. Considering the EOM and going
to the conformal regime, we learn that the Majorana field has dimension
$\Delta_{\chi}=1/q_{\chi}=1/(2q)$. The “spin” correlator $G$ is the square of
the fermion correlator, so in the conformal regime, it will be proportional to
the conformal result with dimension $\Delta=1/q$. This is precisely what we
found numerically for $M=4$ and $q\geq 3$. Note that in the case $q=2$, we do
get SYK-like physics over a wide range of energies at large $M$, but the SYK-
like saddle suffers from an known instability which results in glassy physics
at the lowest energies [13, 25].

###  7.1 Bounds on Glassiness at Large $M$

What happens to the system at large but finite $M$? From our numerics on Spin
SYK in Sec. 6, it seems likely that the model is just SYK-like all the way
down to $M=4$ for sufficiently large $q$. Still, we would like to constrain
possible glassy states.

We can derive a bound on the analog of the Edwards-Anderson order parameter
for $q$-local Gauged Clusters (similar to the recent analysis in [25]). With
the normalization of $G$ chosen above, we have $G(0)=1$. Now introducing
replicas $a=1,\cdots,n$, we have $G_{11}(0)=1$ and we ask how large can
$G_{1b}$ be? This off-diagonal term can be bounded by constraining the
possible expectation values of the $O_{\alpha\beta}$. The best we can do is
choose a state of the fermions in which

$-2i\chi_{1}\chi_{2}=-2i\chi_{3}\chi_{4}=\cdots=-2i\chi_{M-1}\chi_{M}=1,$ (85)

or some permutation thereof. In such a state, $M/2$ of the $\binom{M}{2}$
fermion bilinears have a non-zero expectation value, so

$G_{1b}\leq\frac{1}{M-1}.$ (86)

This formula tells us that glassy physics is suppressed at large $M$. For
$M=4$, this formula reduces to the previous result that $G_{1b}/G_{11}(0)\leq
1/3$.

Suppose that we are at sufficiently large $M$ or $q$ such that $G$ can develop
into its conformal form for some range of $\beta$. $G$ is smallest at
$\tau=\beta/2$ where it is

$G(\beta/2)\sim 1/(\beta J)^{2/q},$ (87)

which only approaches $1/M$ when

$\beta J\sim M^{q/2}.$ (88)

At large $M$ and/or large $q$, this indicates the existence of a
parametrically large window in energy over which the physics can be SYK-like.

###  7.2 Comments on Generic $q$-Local Interactions

We can further generalize the $q$-local Gauged Clusters model by expanding the
set of local operators which can be appear in the $q$-body interaction. So
far, we restricted these operators to be fermion bilinears, but we could allow
more general operators, such as
$\chi_{\alpha}\chi_{\beta}\chi_{\gamma}\chi_{\delta}$. In fact, if we allow
all possible even powers of fermions, we obtain a complete basis for operators
on the $2^{M/2-1}$-dimensional site. We make two comments about this more
generic situation.

First, while it is more complicated to address the case where we add all
possible operators with random couplings of equal variance, it is easier to
discuss adding additional operators perturbatively, i.e. with a random
coupling whose variance is suppressed by a small factor. From the fermion
perspective, this looks like adding higher-$q_{\chi}$ interactions. But such
higher-$q_{\chi}$ interactions are typically irrelevant in the renormalization
group sense. This means we expect the lowest $q_{\chi}$ present in the
Hamiltonian to control the physics at very low energy. Hence, it could be that
adding these terms actually does not strongly modify the very low temperature
physics. In other words, adding generic interactions perturbatively could
preserve the SYK-like physics when it is present at very low energy.

Second, as we add more operators to the Hamiltonian, the corresponding bound
on glassiness improves. For example, suppose we allow every possible single-
site operator in the $q$-body interactions. There are $d^{2}-1$ such operators
(the number generators of SU($d$) and the number of distinct operators built
from even products of fermions), where again $d=2^{M/2-1}$. Let us also
continue to normalize $G$, which now includes a sum over all single site
operators, so that $G(0)=1$. Using the same fermion state as above, we learn
that $d-1$ distinct operators now have a non-vanishing expectation value. This
number is obtained by considering general operators of the form
$(O_{12})^{x_{1}}(O_{34})^{x_{2}}\cdots$, accounting for a double counting
thanks to $(-1)^{F}=1$, and subtracting one for the identity. Hence, we have

$G_{1b}\leq\frac{d-1}{d^{2}-1}=\frac{1}{d+1},$ (89)

which again reduces to $G_{1b}/G_{11}(0)\leq 1/3$ for $d=2$ ($M=4$).

##  8 Physics of the SYK-like Saddle

Here we briefly discuss the low-energy physics when the power-law decaying
$G(\tau)$ controls the quenched free energy. In particular, in SYK the power-
law decay of the Green function is just one part of a whole suite of
interrelated phenomena. Do we obtain the same set of phenomena in Spin SYK?

The general idea is to appeal to the existence of approximate time-
reparameterization symmetry. The first two terms in (16) are manifestly
invariant under reparameterizations, but the third term breaks this symmetry
explicitly. The ground state solution also breaks the symmetry spontaneously,
thus yielding a pattern of symmetry breaking familiar from SYK.

For example, under a general reparameterization $\tau\to f(\tau)$, $\Sigma$ is
modified to

$\Sigma(\tau_{1},\tau_{2})\to
f^{\prime}(\tau_{1})^{1-\Delta}f^{\prime}(\tau_{2})^{1-\Delta}\Sigma(f(\tau_{1}),f(\tau_{2}))$
(90)

with $\Delta=1/q$. Every choice of $h$ now has a reparameterized form,

$h(\tau)\rightarrow f^{\prime}(\tau)^{1-\Delta}h(f(\tau)).$ (91)

If we could rescale the integration variables $s(\tau)$ to
$f^{\prime}(\tau)^{\Delta}s(f(\tau))$, then the $G$ function would transform
as expected. Of course, this is strictly speaking impossible since $s(\tau)$
lives on the sphere and there are additional microscopic terms in the measure.
However, after an appropriate coarse-graining, the coarse-grained $s(\tau)$ is
no longer confined to the sphere. This suggests that the failure of time
reparameterization symmetry is indeed a UV issue.

This means that when we evaluate that action on a reparameterized
configuration,

$f^{\prime}(\tau_{1})^{\Delta}f^{\prime}(\tau_{2})^{\Delta}G(f(\tau_{1}),f(\tau_{2})),$
(92)

the contribution comes from the UV terms which explicitly break the symmetry.
Following the usual analysis, we expect the Schwarzian term to be among these
contributions. If this is the leading term in the action, then the physics is
very close to that of SYK.

Another possibility at finite $q$ and $M$ is that there could be weakly
irrelevant operators which contribute a non-local action that dominates over
the Schwarzian [66]. This is a possibility when the theory has an irrelevant
operator $O$ with dimension $\Delta_{O}\in(1,3/2)$ that appears in the
microscopic action [67]. Whether this possibility is realized depends on the
operator spectrum and on the nature of the UV terms. Even for models where
this possibility is realized, we still have an extensive ground state entropy,
maximal chaos, and the Schwarzian is still present [66].

For example, take $q=3$. This case is glassy for Spin SYK, but for
sufficiently large $M$ in $q$-local Gauged Clusters, the SYK-like solution
should control the quenched free energy over a wide range of temperatures. In
this case, $\Delta=1/3$ and the corresponding long-range fixed point is in the
Gaussian regime. For a Gaussian fixed point, we know the dimensions of all the
operators in the theory and there are seemingly operators with dimensions in
the range $\Delta_{O}\in(1,3/2)$, e.g. $\phi^{4}$ where $\phi$ is some coarse-
grained form of $O_{\alpha\beta}(\tau)$. Nevertheless, it is not clear if
these operators appear in the UV, especially since they do not contain time-
derivatives whereas the coherent state overlaps in the measure are unity
unless the fields vary in time.

These general comments can be elaborated in the various solvable limits we
considered in the paper. Consider first the large $q$ limit. From the
discussion in Section 4 and the numerical solutions of the EOM in Section 5,
we know the SYK-like saddle in this regime hews close to the literal SYK
solution, up to a rescaling of $G$ and a different numerical prefactor
relating $J$ and $\mathcal{J}$. Thus at large $q$ we expect the thermodynamics
to essentially match those of the large $q$ SYK model. This means an SYK-like
extensive zero temperature entropy which approaches the infinite temperature
entropy as $q$ becomes large. We also should have the large $q$ chaos
properties, including maximal chaos at low temperature. Moreover, in the
scaling limit in which $q^{2}/N$ is held fixed as $N\to\infty$, SYK is
solvable [68] and Spin SYK will have the same physics [28, 69].

Similarly, at large $M$ in the $q$-local Gauged Clusters model, we know that
the spin correlation $G$ is determined from $G_{\chi}$, which is exactly the
correlator of a fermionic SYK with $q_{\chi}=2q$. From (25), suitably adjusted
to account for the different normalization of $G$, we learn that

$-\frac{E}{N}\propto\int d\tau(G_{\chi}(\tau))^{2q},$ (93)

which in the conformal regime is identical to the SYK expression up to an
overall factor. The integral is

$-\frac{E}{N}-\text{(UV
part)}\sim\int_{\tau_{c}}^{\beta-\tau_{c}}d\tau\frac{\pi^{2}}{\beta^{2}\sin^{2}\frac{\pi\tau}{\beta}}\sim-\frac{2\pi}{\beta}\cot\frac{\pi\tau_{c}}{\beta}$
(94)

for some UV cutoff $\tau_{c}\sim 1/J$, which yields
$E/N=E_{0}/N+c^{\prime}T^{2}+\cdots$.

##  9 Discussion

In this work, we argued for the existence of SYK-like dynamics in a variety of
bosonic models. Focusing on Spin SYK, we catalogued possible solutions to the
single-replica EOM and then, in lieu of carrying out the minimax procedure, we
used small scale numerical simulations to provide input on the nature of the
quenched saddle. Additional analytical insight was provided by a large $q$
limit. We also introduced and analyzed the $q$-local Gauged Clusters model,
showing that it has SYK-like behavior for sufficiently large $M,q$. Our
numerical results suggest that this behavior persists down to $q=5$ in Spin
SYK. Even when the asymptotically low-temperature physics is glassy, e.g.
$q=2,3,4$, we also gave evidence for an intermediate regime of energies with
SYK-like physics.

We considered a three component model, but there is also interest in a two-
component version as recently discussed in [30]. One can develop an analogous
formulation for the two-component model, with the main difference being that
the effective spin chain has two components (XY) instead of three (XYZ). We
expect the single-replica solutions that we found here will have an analog for
the two-component action. Furthermore, we expect that the two-component model
will also be SYK-like for sufficiently large $q$. We have preliminary results
for the two-component case which we hope to present elsewhere.

These results have significant implications for the structure of quantum field
theories, suggesting that nontrivial CFTs and holographic duals are more
common than we thought. It also has clear implications for near-term quantum
simulators, creating a more viable path to simulating SYK-like physics, and
opening the door to analogous bosonic simulations of more complicated theories
of gravity.

While this paper lays out the basics of $q$-local quantum physics in Spin SYK,
there is much still to do. Some broad questions include:

  1. 1.

Are these results limited to 0+1D systems, or are there other situations where
generic disordered systems have holographic duals? There has been work on
systems with generic local interactions on a graph [70] and in translation
invariant spin chains [71], but much remains to be understood.

  2. 2.

What can be said about fluctuations of the action away from the saddle? Spin
SYK and its cousins have a reparameterization symmetry spontaneously broken
down to $SL(2)$. Does that mean that they have a Schwarzian action that
dominates? Not necessarily, as some systems [66] have non-local effective
actions. It should be possible to understand whether this is realized in Spin
SYK by a careful study of the low-temperature thermodynamics and response
functions.

  3. 3.

For the experimentalists / quantum simulationists, what is the actual two-
point function for Spin SYK at $q=4$ at large system size? In section 5 we
find a replica-diagonal solution consistent with SYK-like physics exists but
the exact diagonalization results suggest a small but non-zero value of
$q_{\textrm{EA}}$. What about for $q$-local Gauged Clusters at $q=3$ for some
modest intermediate $M$, such as $M=6$ or $8$? These questions are all ideal
targets for a quantum simulator.

  4. 4.

For the holographers, since we know that quenched is not equal to annealed in
Spin SYK, what is the gravitational analogue of the annealed solution
(assuming the usual gravity solution is the analogue of the quenched
solution)? How do these solutions interplay with the recent developments
regarding the role of averaging in theories of gravity? Is there is
holographic interpretation of the glassy state for $q=2,3,4$?

  5. 5.

Can we analyze the limit of large local Hilbert space dimension while also
allowing generic local operators in the interaction? If so, is the physics
SYK-like, at least over a wide range of energies, as suggested by the
$q$-local Gauged Clusters model?

Some more model-specific tasks include:

  1. 1.

It would be interesting to carry out a systematic study of the SYK-like
solution at low energies to, for example, definitively establish the form of
the heat capacity and to determine the zero temperature entropy. Perhaps the
argument in [72] can be useful.

  2. 2.

It would likewise be desirable to carry out a full numerical analysis of key
observables in sparse finite $N$ instances, quantities such as the real-time
decay of correlations and out-of-time-order correlations.

  3. 3.

More ambitiously, it would be very interesting to construct replica non-
diagonal saddles and to carry out the minimax procedure. This is work in
progress.

  4. 4.

It would also be very interesting to compute the leading $1/q$ correction [73]
in Spin SYK and the $q$-local Gauged Clusters model as well as the leading
$1/M$ correction in $q$-local Gauged Clusters.

### Acknowledgements

We thank M. Hanada, A. Jevicki, X. Liu, E. Rinaldi, and M. Tezuka for many
discussions. We also thank F. Haehl, C. Baldwin, S. Sachdev, and A. Milekhin
for additional discussions and helpful feedback on the manuscript. We
gratefully acknowledge support from Joint Quantum Institute and from AFOSR
under grant number FA9550-19-1-0360. The numerical data were obtained using
the Zaratan high-performance computing cluster at the University of Maryland,
College Park. We also thank the Perimeter Institute and the It From Qubit
Simons Collaboration for hosting a wonderful Last Hurrah where some of this
work was done.

##  Appendix A The Spin SYK Path Integral

###  A.1 Spin Coherent States

The path integral is derived using spin coherent states. Note that, although
we use the conventional spin language, an instance of the model has no spin
rotation symmetry. For a single spin-1/2 degree of freedom, we take the spin
coherent state to be

$\ket{\Omega}=\begin{bmatrix}\cos\frac{\theta}{2}\\\
\sin\frac{\theta}{2}e^{i\phi}\end{bmatrix}$ (95)

with

$\Omega=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta).$ (96)

The overlap of two spin coherent states is

$\langle\Omega|\Omega^{\prime}\rangle=\cos\frac{\theta}{2}\cos\frac{\theta^{\prime}}{2}+\sin\frac{\theta}{2}\sin\frac{\theta^{\prime}}{2}e^{-i(\phi-\phi^{\prime})}.$
(97)

This set of states yields a resolution the identity ($d\Omega\equiv
d\theta\sin\theta d\phi$):

$\int\frac{d\Omega}{2\pi}|\Omega\rangle\langle\Omega|=\int\frac{d\Omega}{2\pi}\begin{bmatrix}\cos^{2}\frac{\theta}{2}&\cos\frac{\theta}{2}\sin\frac{\theta}{2}e^{-i\phi}\\\
\cos\frac{\theta}{2}\sin\frac{\theta}{2}e^{i\phi}&\sin^{2}\frac{\theta}{2}\end{bmatrix}=I.$
(98)

They can also be used to represent the Pauli matrices,

$\int\frac{d\Omega}{2\pi}|\Omega\rangle\langle\Omega|(3\Omega_{x})=\int\frac{d\Omega}{2\pi}\begin{bmatrix}\cos^{2}\frac{\theta}{2}&\cos\frac{\theta}{2}\sin\frac{\theta}{2}e^{-i\phi}\\\
\cos\frac{\theta}{2}\sin\frac{\theta}{2}e^{i\phi}&\sin^{2}\frac{\theta}{2}\end{bmatrix}(3\sin\theta\cos\phi)=\sigma_{x}.$
(99)
$\int\frac{d\Omega}{2\pi}|\Omega\rangle\langle\Omega|(3\Omega_{y})=\int\frac{d\Omega}{2\pi}\begin{bmatrix}\cos^{2}\frac{\theta}{2}&\cos\frac{\theta}{2}\sin\frac{\theta}{2}e^{-i\phi}\\\
\cos\frac{\theta}{2}\sin\frac{\theta}{2}e^{i\phi}&\sin^{2}\frac{\theta}{2}\end{bmatrix}(3\sin\theta\sin\phi)=\sigma_{y}.$
(100)
$\int\frac{d\Omega}{2\pi}|\Omega\rangle\langle\Omega|(3\Omega_{z})=\int\frac{d\Omega}{2\pi}\begin{bmatrix}\cos^{2}\frac{\theta}{2}&\cos\frac{\theta}{2}\sin\frac{\theta}{2}e^{-i\phi}\\\
\cos\frac{\theta}{2}\sin\frac{\theta}{2}e^{i\phi}&\sin^{2}\frac{\theta}{2}\end{bmatrix}(3\cos\theta)=\sigma_{z}.$
(101)

To demonstrate the machinery, consider a single spin with Hamiltonian
$H=-h\sigma_{x}$. For infinitestimal imaginary time $\Delta\tau$, we write

$e^{-\Delta\tau H}=I+\Delta\tau
h\sigma_{x}=\int\frac{d\Omega}{2\pi}|\Omega\rangle\langle\Omega|(1+\Delta\tau
h(3\Omega_{x})).$ (102)

Diving $\beta$ into $L$ segments of length $\Delta\tau=\beta/L$ and
introducing an $\Omega$ variable for each segment, we have

$\tr(e^{-\beta
H})=\int\left[\prod_{j}\frac{d\Omega_{j}}{2\pi}\langle\Omega_{j+1}|\Omega_{j}\rangle\right]\exp\left(\sum_{j}\Delta\tau
hs_{x}(j)\right)=\int\left[\prod_{j}\frac{d\Omega_{j}}{2\pi}\right]\exp\left(\sum_{j}\Delta(\tau
hs_{x}(j)+i\sin^{2}\theta_{j}\Delta\phi_{j})\right)$ (103)

where $s_{x}(i)=3\Omega_{xi}$. We absorb all the factors in square brackets
into a measure $\mathcal{D}\Omega$ to obtain the compact path integral formula

$\tr(e^{-\beta H})=\int\mathcal{D}\Omega\exp\left(h\int_{0}^{\beta}d\tau
s_{x}(\tau)\right).$ (104)

###  A.2 Replicated Partition Function

Now we give the multi-replica path integral for Spin SYK. For an
infinitestimal imaginary time $\Delta\tau$, we have

$e^{-\Delta\tau H}=I-\Delta\tau\sum_{r_{1}\mu_{1}\cdots
r_{q}\mu_{q}}J_{r_{1}\mu_{1}\cdots
r_{q}\mu_{q}}\sigma_{r_{1}\mu_{1}}\cdots\sigma_{r_{q}\mu_{q}}.$ (105)

Using the coherent state formulas just above and the notation
$s_{\mu}=3\Omega_{\mu}$, we have

$e^{-\Delta\tau
H}=\int\prod_{r}\frac{d\Omega_{r}}{2\pi}|\Omega_{1},\cdots,\Omega_{n}\rangle\langle\Omega_{1},\cdots,\Omega_{n}|\left[1-\Delta\tau\sum_{r_{1}\mu_{1}\cdots
r_{q}\mu_{q}}J_{r_{1}\mu_{1}\cdots r_{q}\mu_{q}}s_{r_{1}\mu_{1}}\cdots
s_{r_{q}\mu_{q}}\right].$ (106)

We have one $\Omega$ variable for each spin and each imaginary time point,
$\Omega_{ri}$. The full measure is

$\mathcal{D}^{N}\Omega=\prod_{r}\left(\prod_{i}\frac{d\Omega_{ri}}{2\pi}\langle\Omega_{ri+1}|\Omega_{ri}\rangle\right).$
(107)

The path integral for the partition function of a single instance is

$Z=\int\mathcal{D}^{N}\Omega\exp\left(-\int d\tau\sum_{r_{1}\mu_{1}\cdots
r_{q}\mu_{q}}J_{r_{1}\mu_{1}\cdots r_{q}\mu_{q}}s_{r_{1}\mu_{1}}(\tau)\cdots
s_{r_{q}\mu_{q}}(\tau)\right).$ (108)

We equivalently write this in terms of a multi-index $I=r_{1}\mu_{1}\cdots
r_{q}\mu_{q}$.

The replicated path integral is obtained from the disorder average of $Z^{m}$,

$\mathbb{E}\left[Z^{m}\right]=\int\mathcal{D}^{N}\Omega_{m}\exp\left(\frac{\mathbb{E}\left[J_{I}^{2}\right]}{2q!}\sum_{I}\left[\sum_{a=1}^{m}\int_{0}^{\beta}d\tau
s^{a}_{r_{1}\mu_{1}}\cdots s^{a}_{r_{q}\mu_{q}}\right]^{2}\right),$ (109)

where now we have included a replica index $a$ on each spin variable (and
modified the measure include $m$ copies, $D\Omega_{m}$). The replica index
runs from $1$ to $m$. The quenched free energy is obtained from the $m\to 0$
limit as discussed in the main text.

We now introduce the $G,\Sigma$ collective fields. The spin correlation
$G_{ab}(\tau_{1},\tau_{2})$ is

$G_{ab}(\tau_{1},\tau_{2})=\frac{1}{N}\sum_{r\mu}s_{r\mu}^{a}(\tau_{1})s_{r\mu}^{b}(\tau_{2}),$
(110)

and combined with the Lagrange multiplier $\Sigma_{ab}$, the path integral
becomes

$\mathbb{E}\left[Z^{m}\right]=\int\mathcal{D}^{N}\Omega_{m}\mathcal{D}\Sigma\mathcal{D}G\exp\left(\frac{J^{2}N}{2q}\sum_{ab}\int
G_{ab}^{q}-\frac{N}{2}\sum_{ab}\int\Sigma_{ab}\left[G_{ab}-\frac{1}{N}\sum_{r\mu}s^{a}_{r\mu}(\tau_{1})s^{b}_{r\mu}(\tau_{2})\right]\right).$
(111)

One important feature is that $G_{a\neq b}(\tau_{1},\tau_{2})$ should be a
constant (possibly zero) since we have separate time-translation invariance in
the individual replicas.

There is another representation of the path integral which decouples the
$s^{a}s^{b}$ term by introducing a fluctuating magnetic field $h_{ra\mu}$ for
each spin in each replica. The fields are independent between spins, and for a
single spin across all replicas, the covariance of $h_{a\mu}(\tau_{1})$ and
$h_{b\nu}(\tau_{2})$ is $\delta_{\mu\nu}\Sigma_{ab}(\tau_{1},\tau_{2})$. This
includes the possibility that $h_{a\mu}(\tau)$ is correlated between replicas.
The path integral is then expressed as

$\mathbb{E}\left[Z^{m}\right]=\int\mathcal{D}^{N}\Omega_{m}\mathcal{D}\Sigma\mathcal{D}G\mathcal{D}^{N}h\exp\left(\frac{J^{2}N}{2q}\sum_{ab}\int
G_{ab}^{q}-\frac{N}{2}\sum_{ab}\int\Sigma_{ab}G_{ab}+\sum_{ra}\int d\tau
h_{ra\mu}(\tau)s^{a}_{r\mu}(\tau)\right),$ (112)

where $\mathcal{D}^{N}h$ is a product of Gaussian measures for $h_{ra\mu}$. If
we now do the $\Omega$ and $h$ integrals, we can write the replicated
partition function as

$\mathbb{E}\left[Z^{m}\right]=\int\mathcal{D}\Sigma\mathcal{D}G\exp\left(\frac{J^{2}N}{2q}\sum_{ab}\int
G_{ab}^{q}-\frac{N}{2}\sum_{ab}\int\Sigma_{ab}G_{ab}+N\log\langle\mathcal{P}e^{\int
d\tau h_{a}\cdot s_{a}}\rangle_{\Sigma_{ab}}\right),$ (113)

The single replica path integral is obtained by setting $m=1$. It is

$\mathbb{E}\left[Z\right]=\int\mathcal{D}^{N}\Omega\mathcal{D}\Sigma\mathcal{D}G\exp\left(\frac{J^{2}N}{2q}\int
G^{q}-\frac{N}{2}\int\Sigma\left[G-\frac{1}{n}\sum_{r\mu}s_{r\mu}(\tau_{1})s_{r\mu}(\tau_{2})\right]\right)$
(114)

or, equivalently,

$\mathbb{E}\left[Z\right]=\int\mathcal{D}^{N}\Omega\mathcal{D}\Sigma\mathcal{D}G\mathcal{D}^{N}h\exp\left(\frac{J^{2}N}{2q}\int
G^{q}-\frac{N}{2}\int\Sigma G+\sum_{r}\int d\tau
h_{r\mu}(\tau)s_{r\mu}(\tau)\right).$ (115)

##  Appendix B Numerical methods

###  B.1 Solving the Equations of Motion

For the data in Sec. 5, we used MATLAB to setup an iterative solver scheme on
a discretized imaginary time circle. The solver has a few components.

First, given a value of $\Sigma$, we need to obtain an estimate of $G$. As
stated in the main text, we use sampling to do this. However, we want to do
importance sampling from the full distribution (single spin partition function
times the Gaussian distribution) of the magnetic fields rather than just
randomly sampling from the Gaussian distribution. To do this, we implemented a
Markov chain Monte Carlo based on the Metropolis-Hastings algorithm. We
generated candidate steps in the magnetic field by sampling from the Gaussian
distribution and multiplying by a small factor, typically of order $.1$ or
$.2$. We then randomly accept or reject the step based on the Metropolis-
Hastings rule. After some number of steps, typically $6$, we use the resulting
fields to generate one sample contributing to $G$. We do not know the auto-
correlation time of our Markov chain, so these parameters are largely ad hoc,
chosen after much experimentation to give reliable solutions.

Second, given a value of $G$, we need to obtain the corresponding $\Sigma$.
This is trivial using the $G$ EOM. We do enforce that $G$ is reflection
symmetric in the above sampling procedure, and we also enforce that $\Sigma$
is symmetric and translation invariant.

Third, we need an update rule. We choose a mix parameter, $x$ and a fraction
$b$, and propose updates via

$G_{\text{prop}}=(1-x)G_{\text{old}}+xG_{\text{new}}$ (116)

where $G_{\text{new}}$ is obtained from the sampling subroutine above. There
is an identical update rule for $\Sigma$. If this update does not increase the
error, then we accept it and continue. If it does increase the error, then we
reduce $x\to bx$ and continue. We always carry out at least one update since
the error is initialized to a large value. Typically, $b=.5$ and $x$ is
initialized to $.5$.

The data in the main text allows for $6$ iterations of this procedure starting
from an initial guess. When not using enough time points or enough samples,
one often finds that no updates after the first are accepted. However, with
sufficient time points and samples, we have found that nearly every update can
be accepted.

The number of iterations and time points are quite modest compared to what is
standard for solutions of the SYK equations of motion. However, the sampling
procedure is currently very expensive. The most expensive runs, e.g. Figure 4,
with $N_{\tau}=300$ time points and $N_{s}=10^{4}$ samples took around $5$
hours to run on a single core. The code was run on the Zaratan cluster at
University of Maryland, College Park.

We are not sure if there is a vastly more efficient method for solving the
equations of motion. However, at least by parallelizing the sampling and
making some other modest improvements, one could probably study significantly
larger $N_{\tau}$ and $N_{s}$.

###  B.2 Finite Size Numerics

For the data in Sec. 6, we used MATLAB and defined instances of the sparse
Spin SYK model as sparse matrices. We then made calls to MATLAB’s built-in
sparse matrix eigenvalue routines to access the lowest few eigenvalues and
eigenvectors of the sparse Hamiltonian. The code was run on the Zaratan
cluster at University of Maryland, College Park.

## References

  * [1] A. J. Bray and M. A. Moore, “Replica theory of quantum spin glasses,” Journal of Physics C: Solid State Physics, vol. 13, p. L655, aug 1980.
  * [2] J. Ye, S. Sachdev, and N. Read, “Solvable spin glass of quantum rotors,” Phys. Rev. Lett., vol. 70, pp. 4011–4014, Jun 1993.
  * [3] S. Sachdev, “Spin glasses enter the quantum regime,” Physics World, vol. 7, p. 25, oct 1994.
  * [4] M. Christos, F. M. Haehl, and S. Sachdev, “Spin liquid to spin glass crossover in the random quantum heisenberg magnet,” Physical Review B, vol. 105, feb 2022.
  * [5] A. Georges, O. Parcollet, and S. Sachdev, “Mean field theory of a quantum heisenberg spin glass,” Physical Review Letters, vol. 85, pp. 840–843, jul 2000.
  * [6] G. Biroli and L. F. Cugliandolo, “Quantum thouless-anderson-palmer equations for glassy systems,” Physical Review B, vol. 64, jun 2001.
  * [7] T. Anous and F. M. Haehl, “The quantum p-spin glass model: a user manual for holographers,” Journal of Statistical Mechanics: Theory and Experiment, vol. 2021, p. 113101, nov 2021.
  * [8] L. Erdős and D. Schröder, “Phase Transition in the Density of States of Quantum Spin Glasses,” Mathematical Physics, Analysis and Geometry, vol. 17, pp. 441–464, Dec. 2014.
  * [9] S. Gharibian, Y. Huang, Z. Landau, and S. W. Shin, “Quantum Hamiltonian Complexity,” arXiv e-prints, p. arXiv:1401.3916, Jan. 2014\.
  * [10] A. W. Harrow and A. Montanaro, “Extremal eigenvalues of local Hamiltonians,” arXiv e-prints, p. arXiv:1507.00739, July 2015.
  * [11] J. P. Keating, N. Linden, and H. J. Wells, “Spectra and Eigenstates of Spin Chain Hamiltonians,” Communications in Mathematical Physics, vol. 338, pp. 81–102, Aug. 2015.
  * [12] Chi-Fang, Chen, A. M. Dalzell, M. Berta, F. G. S. L. Brandão, and J. A. Tropp, “Sparse random Hamiltonians are quantumly easy,” arXiv e-prints, p. arXiv:2302.03394, Feb. 2023.
  * [13] S. Sachdev and J. Ye, “Gapless spin-fluid ground state in a random quantum heisenberg magnet,” Phys. Rev. Lett., vol. 70, pp. 3339–3342, May 1993\.
  * [14] A. Kitaev, “A simple model of quantum holography,” in KITP strings seminar and Entanglement, vol. 12, 2015.
  * [15] J. Maldacena and D. Stanford, “Remarks on the sachdev-ye-kitaev model,” Physical Review D, vol. 94, Nov 2016.
  * [16] V. Rosenhaus, “An introduction to the syk model,” Journal of Physics A: Mathematical and Theoretical, vol. 52, p. 323001, Jul 2019.
  * [17] R. A. Davison, W. Fu, A. Georges, Y. Gu, K. Jensen, and S. Sachdev, “Thermoelectric transport in disordered metals without quasiparticles: The sachdev-ye-kitaev models and holography,” Physical Review B, vol. 95, Apr 2017.
  * [18] A. Gaikwad, L. K. Joshi, G. Mandal, and S. R. Wadia, “Holographic dual to charged syk from 3d gravity and chern-simons,” Journal of High Energy Physics, vol. 2020, Feb 2020.
  * [19] B. Kobrin, Z. Yang, G. D. Kahanamoku-Meyer, C. T. Olund, J. E. Moore, D. Stanford, and N. Y. Yao, “Many-body chaos in the sachdev-ye-kitaev model,” 2020.
  * [20] N. Sorokhaibam, “Phase transition and chaos in charged syk model,” Journal of High Energy Physics, vol. 2020, Jul 2020.
  * [21] W. Fu, D. Gaiotto, J. Maldacena, and S. Sachdev, “Publisher’s note: Supersymmetric sachdev-ye-kitaev models [phys. rev. d , 026009 (2017)],” Physical Review D, vol. 95, mar 2017.
  * [22] F. Sun and J. Ye, “Periodic table of the ordinary and supersymmetric sachdev-ye-kitaev models,” Physical Review Letters, vol. 124, jun 2020\.
  * [23] S. Xu, L. Susskind, Y. Su, and B. Swingle, “A sparse model of quantum holography,” 2020.
  * [24] H. Shackleton, A. Wietek, A. Georges, and S. Sachdev, “Quantum phase transition at nonzero doping in a random $t\text{$-$}j$ model,” Phys. Rev. Lett., vol. 126, p. 136602, Mar 2021.
  * [25] M. Christos, F. M. Haehl, and S. Sachdev, “Spin liquid to spin glass crossover in the random quantum heisenberg magnet,” Phys. Rev. B, vol. 105, p. 085120, Feb 2022.
  * [26] G. Gur-Ari, R. Mahajan, and A. Vaezi, “Does the SYK model have a spin glass phase?,” Journal of High Energy Physics, vol. 2018, p. 70, Nov. 2018\.
  * [27] J. Maldacena, “The Large-N Limit of Superconformal Field Theories and Supergravity,” International Journal of Theoretical Physics, vol. 38, pp. 1113–1133, Jan. 1999.
  * [28] M. Berkooz, P. Narayan, and J. Simón, “Chord diagrams, exact correlators in spin glasses and black hole bulk reconstruction,” Journal of High Energy Physics, vol. 2018, p. 192, Aug. 2018.
  * [29] C. Baldwin and B. Swingle, “Quenched vs annealed: Glassiness from SK to SYK,” Physical Review X, vol. 10, aug 2020.
  * [30] M. Hanada, A. Jevicki, X. Liu, E. Rinaldi, and M. Tezuka, “A model of randomly-coupled Pauli spins,” 9 2023.
  * [31] E. R. Anschuetz, D. Gamarnik, and B. T. Kiani, “Product states optimize quantum $p$-spin models for large $p$,” 2023.
  * [32] L. Foini and J. Kurchan, “Annealed averages in spin and matrix models,” SciPost Physics, vol. 12, mar 2022.
  * [33] T. Tanaka, “Moment problem in replica method,” Interdiscip. Inf. Sci., vol. 13, no. 1, p. 17, 2007.
  * [34] J.-C. Mourrat, “Nonconvex interactions in mean-field spin glasses,” Probab. Math. Phys., vol. 2, no. 2, p. 281, 2021.
  * [35] G. Parisi, “Infinite number of order parameters for spin-glasses,” Phys. Rev. Lett., vol. 43, pp. 1754–1756, Dec 1979.
  * [36] C. L. Baldwin and B. Swingle, “Revisiting the replica trick: Competition between spin glass and conventional order,” Journal of Statistical Physics, vol. 190, jul 2023.
  * [37] K. Usadel and B. Schmitz, “Quantum fluctuations in an ising spin glass with transverse field,” Solid State Communications, vol. 64, no. 6, pp. 975–977, 1987.
  * [38] Y. Y. Goldschmidt, “Solvable model of the quantum spin glass in a transverse field,” Phys. Rev. B, vol. 41, pp. 4858–4861, Mar 1990.
  * [39] A. Vershynina., “Notes on coherent states.”
  * [40] J. SHIBATA and S. TAKAGI, “A NOTE ON (SPIN-) COHERENT-STATE PATH INTEGRAL,” International Journal of Modern Physics B, vol. 13, pp. 107–140, jan 1999.
  * [41] F. J. Dyson, “Existence of a phase-transition in a one-dimensional Ising ferromagnet,” Communications in Mathematical Physics, vol. 12, pp. 91–107, June 1969.
  * [42] M. E. Fisher, S.-k. Ma, and B. G. Nickel, “Critical exponents for long-range interactions,” Phys. Rev. Lett., vol. 29, pp. 917–920, Oct 1972.
  * [43] J. Sak, “Recursion relations and fixed points for ferromagnets with long-range interactions,” Phys. Rev. B, vol. 8, pp. 281–285, Jul 1973.
  * [44] J. Honkonen and M. Y. Nalimov, “Crossover between field theories with short-range and long-range exchange or correlations,” Journal of Physics A: Mathematical and General, vol. 22, p. 751, mar 1989.
  * [45] M. C. Angelini, G. Parisi, and F. Ricci-Tersenghi, “Relations between short-range and long-range ising models,” Physical Review E, vol. 89, jun 2014.
  * [46] M. F. Paulos, S. Rychkov, B. C. van Rees, and B. Zan, “Conformal invariance in the long-range ising model,” Nuclear Physics B, vol. 902, pp. 246–291, jan 2016.
  * [47] C. Behan, L. Rastelli, S. Rychkov, and B. Zan, “Long-range critical exponents near the short-range crossover,” Physical Review Letters, vol. 118, jun 2017.
  * [48] G. Tarnopolsky, “Large $q$ expansion in the sachdev-ye-kitaev model,” Phys. Rev. D, vol. 99, p. 026010, Jan 2019.
  * [49] D. R. Grempel and M. J. Rozenberg, “Fluctuations in a quantum random heisenberg paramagnet,” Phys. Rev. Lett., vol. 80, pp. 389–392, Jan 1998\.
  * [50] A. M. García-García, Y. Jia, D. Rosa, and J. J. M. Verbaarschot, “Sparse sachdev-ye-kitaev model, quantum chaos, and gravity duals,” Phys. Rev. D, vol. 103, p. 106002, May 2021.
  * [51] M. Tezuka, O. Oktay, E. Rinaldi, M. Hanada, and F. Nori, “Binary-coupling sparse sachdev-ye-kitaev model: An improved model of quantum chaos and holography,” Physical Review B, vol. 107, feb 2023.
  * [52] G. Parisi, “Order parameter for spin-glasses,” Phys. Rev. Lett., vol. 50, pp. 1946–1948, Jun 1983.
  * [53] H. Sompolinsky and A. Zippelius, “Relaxational dynamics of the edwards-anderson model and the mean-field theory of spin-glasses,” Phys. Rev. B, vol. 25, pp. 6860–6875, Jun 1982.
  * [54] H.-J. Sommers, “What is the static edwards-anderson order parameter in the mean-field theory of spin glasses?,” Journal de Physique Lettres, vol. 43, no. 21, pp. 719–725, 1982.
  * [55] M. Winer, R. Barney, C. L. Baldwin, V. Galitski, and B. Swingle, “Spectral form factor of a quantum spin glass,” Journal of High Energy Physics, vol. 2022, sep 2022.
  * [56] R. Barney, M. Winer, C. L. Baldwin, B. Swingle, and V. Galitski, “Spectral statistics of a minimal quantum glass model,” SciPost Physics, vol. 15, sep 2023.
  * [57] F. Haake, Quantum Signatures of Chaos.  Springer Series in Synergetics, Springer Berlin Heidelberg, 2010.
  * [58] M. V. Berry and M. Tabor, “Level clustering in the regular spectrum,” Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, vol. 356, no. 1686, pp. 375–394, 1977.
  * [59] F. J. Dyson, “Statistical theory of the energy levels of complex systems. iii,” Journal of Mathematical Physics, vol. 3, no. 1, pp. 166–175, 1962\.
  * [60] E. Wigner and J. Griffin, Group Theory and Its Application to the Quantum Mechanics of Atomic Spectra.  Pure and applied Physics, Academic Press, 1959.
  * [61] M. Mehta, Random Matrices.  ISSN, Elsevier Science, 2004.
  * [62] O. Bohigas, M. J. Giannoni, and C. Schmit, “Characterization of chaotic quantum spectra and universality of level fluctuation laws,” Phys. Rev. Lett., vol. 52, pp. 1–4, Jan 1984.
  * [63] A. M. García-García and J. J. M. Verbaarschot, “Analytical spectral density of the Sachdev-Ye-Kitaev model at finite N,” Physical Review D, vol. 96, p. 066012, Sept. 2017.
  * [64] J. S. Cotler, G. Gur-Ari, M. Hanada, J. Polchinski, P. Saad, S. H. Shenker, D. Stanford, A. Streicher, and M. Tezuka, “Black holes and random matrices,” Journal of High Energy Physics, vol. 2017, p. 118, May 2017.
  * [65] G. Bentsen, Y. Gu, and A. Lucas, “Fast scrambling on sparse graphs,” Proceedings of the National Academy of Science, vol. 116, pp. 6689–6694, Apr. 2019.
  * [66] A. Milekhin, “Coupled sachdev-ye-kitaev models without schwartzian dominance,” 2021.
  * [67] J. Maldacena, D. Stanford, and Z. Yang, “Conformal symmetry and its breaking in two dimensional Nearly Anti-de-Sitter space,” arXiv e-prints, p. arXiv:1606.01857, June 2016.
  * [68] M. Berkooz, M. Isachenkov, V. Narovlansky, and G. Torrents, “Towards a full solution of the large N double-scaled SYK model,” Journal of High Energy Physics, vol. 2019, p. 79, Mar. 2019.
  * [69] M. Berkooz, N. Brukner, S. F. Ross, and M. Watanabe, “Going beyond ER=EPR in the SYK model,” Journal of High Energy Physics, vol. 2022, p. 51, Aug. 2022.
  * [70] R. Movassagh, “Generic local hamiltonians are gapless,” Phys. Rev. Lett., vol. 119, p. 220504, Nov 2017.
  * [71] M. Lemm, “Gaplessness is not generic for translation-invariant spin chains,” Phys. Rev. B, vol. 100, p. 035113, Jul 2019.
  * [72] A. Milekhin, “On minimal residual entropy in non-Fermi liquids,” arXiv e-prints, p. arXiv:2207.01588, July 2022.
  * [73] G. Tarnopolsky, “On large $q$ expansion in the Sachdev-Ye-Kitaev model,” arXiv e-prints, p. arXiv:1801.06871, Jan. 2018.