aainstitutetext: Department of Physics and Astronomy, Stony Brook University,
Stony Brook, New York 11794, USA

# Chaos on the hypercube

Yiyang Jia a and Jacobus J. M. Verbaarschot<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We analyze the spectral properties of a $d$-dimensional HyperCubic (HC)
lattice model originally introduced by Parisi. The U(1) gauge links of this
model give rise to a magnetic flux of constant magnitude $\phi$ but random
orientation through the faces of the hypercube. The HC model, which also can
be written as a model of $2d$ interacting Majorana fermions, has a spectral
flow that is reminiscent of Maldacena-Qi (MQ) model, and its spectrum at
$\phi=0$, actually coincides with the coupling term of the MQ model. As was
already shown by Parisi, at leading order in $1/d$ , the spectral density of
this model is given by the density function of the Q-Hermite polynomials,
which is also the spectral density of the double-scaled Sachdev-Ye-Kitaev
model. Parisi demonstrated this by mapping the moments of the HC model to
$Q$-weighted sums on chord diagrams. We point out that the subleading moments
of the HC model can also be mapped to weighted sums on chord diagrams, in a
manner that descends from the leading moments. The HC model has a magnetic
inversion symmetry that depends on both the magnitude and the orientation of
the magnetic flux through the faces of the hypercube. The spectrum for fixed
quantum number of this symmetry exhibits a transition from regular spectra at
$\phi=0$ to chaotic spectra with spectral statistics given by the Gaussian
Unitary Ensembles (GUE) for larger values of $\phi$. For small magnetic flux,
the ground state is gapped and is close to a Thermofield Double (TFD) state.

##  1 Introduction

Many-body chaos has attracted a great deal of attention in recent years. In
particular, the study of the Sachdev-Kitaev-Ye (SYK) model sachdev1993 ;
Kitaev2015  formerly known as the two-body random ensemble mon1975 , has
greatly improved our understanding of the relationship between many-body
chaos, disorder and spectral properties of the underlying Hamiltonian (see
brody1981 ; benet2003 ; Borgonovi:2016a ; Borgonovi:2019mrk  for reviews and
recent work). One of the main conclusions is that the SYK model is a non-Fermi
liquid with a level density that increases exponentially with the volume
rather than a power of the volume for a Fermi liquid. A direct consequence is
that the zero temperature limit of this model has an nonzero extensive entropy
Georges_2001 . For the same reason, the SYK model can be used to address
questions related to understanding micro-states and entropy of black holes
sachdev2015 .

There are different ways to measure the chaotic properties of the SYK model.
The short-time behavior of the Out-of-Time-Order Correlator (OTOC), which in
the classical limit describes the exponential divergence of classical
trajectories, was shown Maldacena2016  to saturate the chaos bound
Maldacena2015 . This is also expected to be the case for black holes, and was
one of the main reasons for the excitement for the SYK model. The paradigm of
quantum chaos, though, is that spectral correlations are given by Random
Matrix Theory (RMT), which is known as the Bohigas-Giannoni-Schmidt conjecture
bohigas1984 ; Seligman1984 . Indeed this was confirmed by numerical and
analytical studies of the SYK model Garcia-Garcia:2016mno ; Cotler2016 ;
Saad:2018bqo ; Altland:2017eao ; Jia:2019orl . One issue that has come forth
in the study of the SYK model is to what extent the disorder contributes to
its chaotic properties. It has been known for a long time Flores_2001  that
level fluctuations at the scale of many level spacings are dominated by
fluctuations of the width of the spectrum going from one disorder realization
to the next. In the time domain, these fluctuations Altland:2017eao ; Garcia-
Garcia:2018ruf ; Saad:2018bqo ; Gharibyan:2018jrp ; Jia:2019orl  give rise to
a peak at very short times in the connected spectral form factor. This peak
should not be confused with the peak due to the disconnected part of the
spectral form factor which is many orders of magnitude larger. Fluctuations of
other low-order moments also give significant contributions to long range
spectral fluctuations. For an SYK system of $N$ Majorana fermions, the
deviation from Random Matrix Theory are described by the covariance matrix of
the first $O(N)$ moments.111We do not claim certainty on the estimate $O(N)$,
since it is inferred by observing limited numerics. In fact in Jia:2019orl  by
the present authors, another estimate of $O(N\log N)$ was derived
analytically, but that estimate was also based on a crude estimate. This gives
an estimate of $2^{N/2}/N$ for the spectral range of RMT fluctuations or a
times scale of $N2^{-N/2}$ beyond which the spectral form factor is given by
RMT.

It has been argued that the disorder is not important for the correlation
functions and thermodynamics of the SYK model Witten:2016iux  which also has
been confirmed by melonic models which have similar properties in the absence
of disorder Klebanov:2016xxf ; Klebanov:2019jup ; Kim:2019upg ;
Krishnan:2017lra . In this paper we study an SYK-like model with much less
disorder than the SYK model. This is the hypercube $U(1)$ lattice model in $d$
dimensions originally introduced by Parisi Parisi:1994jg ; Marinari:1995jwr ;
Capelli:1997pm ; Colomo:2001a ; Colomo:2002rk  as a model for an array of
Josephson junctions. This model has a magnetic flux of constant magnitude
through each of the faces of a $d$-dimensional hypercube, and only the sign of
the flux through each face is random. In spite of the $U(1)$ disorder on the
links, the first six moments of the spectral density do not depend on the
disorder realization, and the scale fluctuations that limit the agreement with
random matrix theory are absent in this model. Experience with the 2+4-body
SYK model shows Garcia-Garcia:2017bkg ; Nosaka:2018iat ; Nosaka:2019tcx  that
although the two-body term is relevant, the model still remains chaotic, and
also in the hypercubic model we expect to find spectra correlated according to
Random Matrix Theory.

In this paper we study the chaotic properties of Parisi’s Hypercubic (HC)
model. The magnitude of the flux or the Wilson loop is is parameterized by
$\phi$. At $\phi=0$ the Hamiltonian is given by the adjacency matrix of the
hypercube graph which is integrable and coincides with the coupling term of
the Maldacena-Qi model Maldacena:2018lmt . The spectral flow as a function of
$\phi$ is also similar to that of the Maldacena-Qi model, and exhibits an
integrable-to-chaos transition. In addition, the hypercubic model has a
previously unknown discrete symmetry, which is a variant of the magnetic
translation symmetry Zak:1964 ; Rammal:1990 ; Wiegmann:1994js , and is
reminiscent of the discrete symmetry of the Maldacena-Qi model. Understanding
of the exact symmetries is essential for a statistical analysis of the
spectral correlations of this model.

As was already noted by Parisi, the average spectral density for large $d$ is
well approximated by the $Q$-Hermite spectral density with $Q=\cos\phi$. This
also is the case for the double scaling limit with $q^{2}/N$ fixed
($Q=e^{-2q^{2}/N}$) for the $q$-body SYK model of $N$ interacting Majorana
fermions Bagrets2016 ; Garcia-Garcia:2017pzl ; Cotler2016 ; Jia:2018ccl ;
Berkooz:2018qkz ; Berkooz:2018jqr . For $\phi>\pi/2$, $Q$ becomes negative and
spectrum splits into two bands, which also happens for the supercharge of the
supersymmetric SYK model Fu:2016vas ; Garcia-Garcia:2018ruf ; Kanazawa:2017dpd
. The spectral fluctuations of the HC model from one realization to the next
are quite different from those of the SYK model. In the SYK model these
fluctuations result from the covariance of the first $O(N)$ moments, decouple
from the RMT fluctuations quite well, and can be eliminated Flores_2001 ;
Altland:2017eao ; Garcia-Garcia:2018ruf ; Saad:2018bqo ; Gharibyan:2018jrp ;
Jia:2019orl . For HC model, which can also formulated in terms of gamma
matrices in $2d$ dimensions, the fluctuations due to the first few moments are
suppressed, but higher moments seem to contribute significantly to the
deviation from RMT level statistics. The scale of these fluctuations does not
seem to separate well from the scale of the RMT fluctuations.

This model has a gap that seems to remain in the thermodynamical limit for
$\phi<\pi/2$. Therefore the ground state entropy vanishes at zero temperature.
Since for zero flux the model coincides with the coupling Hamiltonian of the
Maldacena-Qi model, the ground state is also given by a ThermoField Double
(TFD) state. However, contrary to the Maldacena-Qi model, the overlap with the
TFD state decreases considerably for nonzero magnetic flux.

This paper is organized as follows. In section 2 we introduce Parisi’s
hypercubic model which, as is explained in section 3, can also be expressed as
a sum of tensor products of Pauli matrices. The novel discrete symmetry of
this model is discussed in section 4. In section 5 we show that the first six
moments of this model do not depend on the disorder realization. Numerical
results for the spectral density and spectral correlations are presented in
section 6. Both the number variance and the spectral form factor are compared
to random matrix results. The ground state wave function is compared to the
TFD state in section 7 and concluding remarks are made in section 8. Several
technical results are worked out in two appendices. In appendix A we calculate
the fourth and sixth moments of the Hamiltonian in a tensor product
representation, respectively. The connection between chord diagrams and the
leading large $d$ moments of the Hamiltonian is explained in appendix B, where
we also demonstrate how subleading moments arise from chord diagram
considerations.

##  2 Parisi’s hypercubic model

Parisi Parisi:1994jg  studied a disordered $U(1)$ lattice gauge model on a
$d$-dimensional Euclidean hypercube. The sites of this model are represented
by $d$-dimensional vectors $\vec{x}$ with components $x_{\mu}\in\\{0,1\\}$.
The model considers a constant magnetic field such that the fluxes through all
faces of the hypercube have the same magnitude $\phi$, but with random
orientations. That is, we have the field strength tensor

$F_{\mu\nu}=\phi S_{\mu\nu},$ (1)

where $S_{\mu\nu}$ is an antisymmetric tensor with random entries $\pm 1$ with
equal probabilities. Hence we are dealing with a finite ensemble with
$2^{\binom{d}{2}}$ disorder realizations. We can work in the axial gauge so
that the link variables are given by

$U_{\mu}(\vec{x})=e^{i\phi\sum_{\nu=1}^{\mu-1}S_{\mu\nu}x_{\nu}},$ (2)

which is the phase we associate with the link emanating from site $\vec{x}$
along the $\mu$-th direction. Note the sum is over all the $\nu$’s with
$\nu<\mu$, and if $\mu=1$ we define $U_{1}(\vec{x})=1$. We wish to study a
Hamiltonian describing a particle hopping on the lattice sites through the
lattice links, and picking up a phase of the corresponding link variable. In
terms of matrix elements, the Hamiltonian $H$ has the form

$H_{\vec{x},\vec{y}}=\sum_{\mu}\left[\delta_{\vec{x}+\hat{e}^{\mu},\vec{y}}\
U_{\mu}(\vec{x})+\delta_{\vec{x}-\hat{e}^{\mu},\vec{y}}\
U^{*}_{\mu}(\vec{x})\right],$ (3)

where $\hat{e}^{\mu}$ is the unit basis vector in the $\mu$-th direction. When
$\phi=0$, this Hamiltonian becomes the adjacency matrix of the hypercube as a
graph. We remark that Parisi was originally interested in the second quantized
Hamiltonian

$\sum_{\vec{x},\vec{y}}\varphi_{\vec{x}}\ H_{\vec{x},\vec{y}}\
\varphi_{\vec{y}},$ (4)

where $\varphi$ is a scalar quantum field. However, in this paper we take a
first quantized view and concern ourselves with the $H$ defined in equation
(3), and the wave functions live in $\mathbb{C}^{2^{d}}$.

Let us be very explicit on how to write the Hamiltonian matrix as a two-
dimensional array of numbers: since $\vec{x}=(x_{1},x_{2},\ldots,x_{d})$ is a
string of $0$ and $1$’s of length $d$, we can naturally think of $\vec{x}$ as
the binary representation of some integer between $0$ and $2^{d}-1$. Shifting
this correspondence by one, we can represent any integer
$m\in\\{1,2,3,\ldots,2^{d}\\}$ through the relation

$[m-1]_{2}=\overline{x_{d}x_{d-1}\ldots x_{1}},$ (5)

where $[m-1]_{2}$ denotes the number $m-1$ in the binary representation, and
$\overline{x_{d}x_{d-1}\ldots x_{1}}$ denotes $x_{d}x_{d-1}\ldots x_{1}$ as a
string of digits. We will use $m$ to index the matrix entries. Note we use the
reverse order of $(x_{1},x_{2},\ldots,x_{d})$ to represent binary digits
because we wish contributions from lower dimensions to appear as the upper-
left of the matrix. For example, with these conventions we have

$\begin{split}H_{d=1}&=\begin{pmatrix}0&\ 1\\\ 1&\ 0\end{pmatrix},\\\
H_{d=2}&=\begin{pmatrix}0&1&1&0\\\ 1&0&0&\ e^{i\phi S_{21}}\\\ 1&0&0&1\\\ 0&\
e^{-i\phi S_{21}}&1&0\end{pmatrix},\end{split}$ (6)

and so on.

The Hamiltonian can be obtained recursively:

$\displaystyle H_{d}=\left(\begin{array}[]{cc}H_{d-1}&C_{d-1}\\\
(C_{d-1})^{-1}&H_{d-1}\end{array}\right),$ (9)

where $C_{d-1}$ is a diagonal unitary matrix with entries

$(C_{d-1})_{k,k}=e^{i\phi\sum_{\nu=1}^{d-1}S_{d\nu}x_{\nu}(k)},$ (10)

where $x_{\nu}(k)$ is the $\nu$-th digit of $[k-1]_{2}$, as defined in
equation (5). We can verify that the following relation holds:

$\displaystyle\left(C_{d-1}\right)_{k,k}\left(C_{d-1}\right)_{2^{d-1}+1-k,\;2^{d-1}+1-k}=e^{i\phi\sum_{\nu=1}^{d-1}S_{d\nu}},$
(11)

which will be useful for section 4. For later convenience, we also introduce
the notation

$\displaystyle S_{\rho}:=\sum_{\nu=1}^{\rho-1}S_{\rho\nu},\quad S_{1}:=0,$
(12)

so that the right-hand side of equation (11) is simply $e^{i\phi S_{d}}$.

##  3 Tensor product representation of the Hamiltonian

Since the interaction between two sites can be written in terms of the Pauli
matrix $\sigma_{1}$, it is not surprising that the Hamiltonian can be
expressed in terms of tensor products of Pauli $\sigma$ matrices. For $d=2$ it
is clear from equation (6) that

$\displaystyle H_{2}$ $\displaystyle=$
$\displaystyle\sigma_{0}\otimes\sigma_{1}+\sigma_{2,(0,0)}\otimes\rho_{0}+\sigma_{2,(1,0)}\otimes\rho_{1}$

with

$\displaystyle\rho_{0}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}1&0\\\
0&0\end{array}\right),\qquad\rho_{1}=\left(\begin{array}[]{cc}0&0\\\
0&1\end{array}\right),\qquad\sigma_{0}=\left(\begin{array}[]{cc}1&0\\\
0&1\end{array}\right),$ (20) $\displaystyle\sigma_{1}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}0&1\\\
1&0\end{array}\right),\qquad\sigma_{d,\vec{x}}=\left(\begin{array}[]{cc}0&e^{i\phi\sum_{\nu=1}^{d-1}x_{\nu}S_{d\nu}}\\\
e^{-i\phi\sum_{\nu=1}^{d-1}x_{\nu}S_{d\nu}}&0\end{array}\right),$ (25)

where $\vec{x}=(x_{1},\ldots,x_{d})$. Notice that the definition of
$\sigma_{d,\vec{x}}$ does not depend on the last component $x_{d}$ of
$\vec{x}$, for example we have

$\sigma_{2,(0,0)}=\sigma_{2,(0,1)}=\begin{pmatrix}0&1\\\
1&0\end{pmatrix}=\sigma_{1}.$ (26)

For higher dimensions we have

$\displaystyle H_{3}$ $\displaystyle=$ $\displaystyle\sigma_{0}\otimes
H_{2}+\sum_{x_{1},x_{2}}\sigma_{3,\vec{x}}\otimes\rho_{x_{1}}\otimes\rho_{x_{2}},$
$\displaystyle H_{4}$ $\displaystyle=$ $\displaystyle\sigma_{0}\otimes
H_{3}+\sum_{x_{1},x_{2},x_{3}}\sigma_{4,\vec{x}}\otimes\rho_{x_{1}}\otimes\rho_{x_{2}}\otimes\rho_{x_{3}},$
(27)

and in general we have

$\displaystyle H_{d}=\sigma_{0}\otimes
H_{d-1}+\sum_{x_{1},x_{2},\cdots,x_{d-1}}\sigma_{d,\vec{x}}\otimes\rho_{x_{1}}\otimes\rho_{x_{2}}\otimes\cdots\otimes\rho_{x_{d-1}},$
(28)

where

$\sum_{x_{1},x_{2},\cdots,x_{d-1}}:=\sum_{x_{1}=0}^{1}\sum_{x_{2}=0}^{1}\cdots\sum_{x_{d-1}=0}^{1}.$
(29)

###  3.1 Hamiltonian as system of interacting Majorana fermions

Since the Hamiltonian is a sum of tensor products of Pauli matrices, it is
natural to express the Hamiltonian as a sum of products of $\gamma$ matrices,
which then can be interpreted and the Hamiltonian for a system of $2d$
Majorana fermions. The simples case is $\phi=0$. Then the Hamiltonian is just
the adjacency matrix of a hypercube graph. In the tensor product
representation it is given by

$\displaystyle H_{d}(\phi=0)=\sigma_{0}\otimes
H_{d-1}(\phi=0)+\sigma_{1}\otimes\overbrace{\sigma_{0}\otimes\cdots\otimes\sigma_{0}}^{d-1},\quad
H_{1}(\phi=0)=\sigma_{1}.$ (30)

If we introduce the gamma matrices

$\displaystyle\gamma_{k}^{L}$ $\displaystyle=$
$\displaystyle\overbrace{\sigma_{1}\otimes\cdots\otimes\sigma_{1}}^{k-1}\otimes\sigma_{3}\otimes\overbrace{\sigma_{0}\otimes\cdots\otimes\sigma_{0}}^{d-k},$
$\displaystyle\gamma_{k}^{R}$ $\displaystyle=$
$\displaystyle\overbrace{\sigma_{1}\otimes\cdots\otimes\sigma_{1}}^{k-1}\otimes\sigma_{2}\otimes\overbrace{\sigma_{0}\otimes\cdots\otimes\sigma_{0}}^{d-k},$
(31)

then the Hamiltonian can be written as

$\displaystyle H_{d}(\phi=0)=i\sum_{k=1}^{d}\gamma_{k}^{L}\gamma_{k}^{R}.$
(32)

This is exactly the interaction term in the Maldacena-Qi model
Maldacena:2018lmt . This interaction term was shown Garcia-Garcia:2019poj  to
have the spectrum

$-d+2i,\ i=0,1,\ldots,d,$ (33)

with degeneracies

$\binom{d}{i},\ i=0,1,\ldots,d.$ (34)

Indeed this is also the well-known spectrum for the hypercube adjacency
matrix. At $\phi\neq 0$ most other terms contributing to the Hamiltonian
couple the L and R spaces, which makes this model quite different from the
Maldacena-Qi model. In addition, interaction terms among any number of
$\gamma$ matrices appear in the Hamiltonian, which makes the Hamiltonian look
very complicated.

##  4 Symmetries

###  4.1 Sublattice symmetry

Since the hypercube is a bi-partite lattice, the lattice links only connect
one sublattice to the other, we conclude that the Hamiltonian (3) has a
sublattice symmetry

$\displaystyle\\{\Gamma_{5},H\\}=0.$ (35)

so that all eigenvalues appear in pairs $\pm\lambda_{k}$. In the tensor
product representation described in section 3, $\Gamma_{5}$ has the simple
form of a tensor product of $\sigma_{3}$’s. This is because each term
contributing to $H_{d}$ in equation (28) contains exactly one off-diagonal
$\sigma$ matrix in the tensor product. Therefore,

$\displaystyle\\{\overbrace{\sigma_{3}\otimes\cdots\otimes\sigma_{3}}^{d},H_{d}\\}=0,$
(36)

which proves the sublattice symmetry of the Hamiltonian.

###  4.2 Magnetic inversion symmetry

Since the field strength is constant in space and is a two-form, it is
invariant under inversion

$\vec{x}\to\vec{x}^{\,c}:=(1-x_{1},1-x_{2},\ldots,1-x_{d}).$ (37)

We choose the inverted coordinates to be $1-x_{k}$ instead of $-x_{k}$ so that
the hypercube remains invariant too. Therefore, we expect a symmetry $A_{d}$
of the system acting on wave functions as

$A_{d}\psi(\vec{x})=\Omega_{d}(\vec{x})\psi\left(\vec{x}^{\,c}\right),$ (38)

and we seek a position-dependent phase factor $\Omega_{d}(\vec{x})$ so that
$[A_{d},H_{d}]=0$. Its global phase is still ambiguous, which can be fixed by
requiring $A_{d}^{2}=\mathbb{1}_{2^{d}\times 2^{d}}$ as a phase convention. We
claim that for $A_{d}$ to commute with the Hamiltonian, we have to make the
choice

$\Omega_{d}(x)=\exp\left(i\frac{\phi}{2}\sum_{\rho=2}^{d}S_{\rho}\right)\exp\left({i\phi\sum_{\nu>\mu}S_{\mu\nu}x_{\nu}}\right),$
(39)

where $S_{\rho}$ is defined as in equation (12) and

$\sum_{\nu>\mu}:=\sum_{\nu=2}^{d}\sum_{\mu=1}^{\nu-1}.$ (40)

Now equations (37) and (39) fix $A_{d}$ unambiguously. We can write $A_{d}$
explicitly as a matrix through the recursion relation:

$\displaystyle A_{d}=\left(\begin{array}[]{cc}0&e^{\frac{i}{2}\phi
S_{d}}A_{d-1}\\\ e^{-\frac{i}{2}\phi S_{d}}A_{d-1}&0\end{array}\right),\qquad
A_{1}=\sigma_{1}.$ (43)

Note that $A_{d}$ is a Hermitian anti-diagonal matrix. By induction we easily
check that $A_{d}^{2}=\mathbb{1}$ so that its eigenvalues can only be $\pm 1$.
We will call $A_{d}$ the magnetic inversion, because the operator implements a
spatial inversion and is a function of the magnetic field. Let us remark that
although we only wanted to implement a inversion in space, since we wrote down
the Hamiltonian in a specific gauge (in our case the axial gauge along $x_{1}$
direction), spatial transformations may not always respect the gauge
condition. The position-dependent phase factor is the price to pay to stay in
the same gauge. This is analogous to the more familiar case of magnetic
translations. In fact $A_{d}$ can be viewed as a magnetic translation if we
view the inversion (37) as a translation mod 2 along the longest body diagonal
of the hypercube:

$1-x_{k}=(1+x_{k})\text{ mod 2},\text{ for all }k=1,2,\ldots,d.$ (44)

Then the position-dependent phase factor in equation (39) is exactly the
standard phase factor for the corresponding magnetic translation
Sekiguchi:2008kw .

We now prove $A_{d}$ is indeed a symmetry by induction. The commutator
$[A_{d},H_{d}]$ is given by

$\displaystyle\begin{pmatrix}e^{\frac{i}{2}\phi
S_{d}}A_{d-1}(C_{d-1})^{-1}-e^{-\frac{i}{2}\phi
S_{d}}C_{d-1}A_{d-1}&e^{\frac{i}{2}\phi S_{d}}[A_{d-1},H_{d-1}]\\\
e^{-\frac{i}{2}\phi S_{d}}[A_{d-1},H_{d-1}]&e^{-\frac{i}{2}\phi
S_{d}}A_{d-1}C_{d-1}-e^{\frac{i}{2}\phi
S_{d}}(C_{d-1})^{-1}A_{n-1}\end{pmatrix}.\quad$ (45)

By induction hypothesis $[A_{d-1},H_{d-1}]=0$ which is satisfied for $d=2$
because $H_{1}=A_{1}=\sigma_{1}$, so we only have to worry about the diagonal
blocks. We remind the readers that $C_{d-1}$ is diagonal and $A_{d-1}$ is
anti-diagonal, so their product is anti-diagonal. So let us look at the only
matrix elements that are possibly nonzero:

$\displaystyle\begin{split}&\left[e^{\frac{i}{2}\phi
S_{d}}A_{d-1}(C_{d-1})^{-1}\right]_{k,2^{d-1}+1-k}-\left[e^{-\frac{i}{2}\phi
S_{d}}C_{d-1}A_{d-1}\right]_{k,2^{d-1}+1-k}\\\ =&e^{\frac{i}{2}\phi
S_{d}}\left[A_{d-1}\right]_{k,2^{d-1}+1-k}\left[C_{d-1}\right]_{2^{d-1}+1-k,2^{d-1}+1-k}^{-1}-e^{-\frac{i}{2}\phi
S_{d}}\left[C_{d-1}\right]_{kk}\left[A_{d-1}\right]_{k,2^{d-1}+1-k}\\\
=&e^{-\frac{i}{2}\phi
S_{d}}\left[A_{d-1}\right]_{k,2^{d-1}+1-k}\left[C_{d-1}\right]_{2^{d-1}+1-k,2^{d-1}+1-k}^{-1}\left(e^{i\phi
S_{d}}-\left[C_{d-1}\right]_{kk}\left[C_{d-1}\right]_{2^{d-1}+1-k,2^{d-1}+1-k}\right)\\\
=&0,\end{split}$

where for the last equality we used equation (11) and this completes the
proof.

Since the symmetry operator is an anti-diagonal matrix, an orthogonal set of
eigenvectors is given by $(0,\cdots,0,b_{k},0,\dots,0,\pm
b_{2^{d}+1-k},0,\cdots,0)$, where the $b_{k}$ are the anti-diagonal matrix
elements. For the symmetry operator $A_{d}$ we have that
$b_{2^{d}+1-k}=b_{k}^{*}$. These eigenvectors can be used to construct the
unitary matrix that brings the Hamiltonian into a block-diagonal form where
the two blocks correspond to the $\pm 1$ eigenvalues of $A_{d}$. For our
numerical results to be discussed below, we block-diagonalize the Hamiltonian
this way.

We can also discuss the magnetic inversion symmetry in the tensor product
representation. If we define the unitary Hermitian $2\times 2$ matrix

$\displaystyle K_{d}:=\left(\begin{array}[]{cc}0&e^{i\frac{\phi}{2}S_{d}}\\\
e^{-i\frac{\phi}{2}S_{d}}&0\end{array}\right),\quad
K_{1}:=\sigma_{1}=\begin{pmatrix}0&1\\\ 1&0\end{pmatrix},$ (48)

then it is clear from equation (43) that the magnetic inversion can be written
as

$A_{d}=K_{d}\otimes K_{d-1}\otimes\cdots\otimes K_{1}.$ (49)

We can check that

$\displaystyle K_{d}\sigma_{d,\vec{x}}K_{d}^{-1}$ $\displaystyle=$
$\displaystyle\sigma_{d,\vec{x}^{c}},$ $\displaystyle
K_{d}\rho_{x_{k}}K_{d}^{-1}$ $\displaystyle=$ $\displaystyle\rho_{x_{k}^{c}},$
(50)

where $x_{k}^{c}=1-x_{k}$ and
$\vec{x}^{\,c}=(x_{1}^{c},x_{2}^{c}\ldots,x_{d}^{c})$, as defined in equation
(37). Now we can prove $A_{d}H_{d}A_{d}^{-1}=H_{d}$ from induction again: the
inductive hypothesis takes care of the first term on the right-hand side of
equation (28), and the second term becomes

$\sum_{x_{1},x_{2},\cdots,x_{d-1}}\sigma_{d,\vec{x}^{c}}\otimes\rho_{x_{1}^{c}}\otimes\rho_{x_{2}^{c}}\otimes\cdots\otimes\rho_{x_{d-1}^{c}},$
(51)

but we can freely re-index the summation as

$\sum_{x_{1}^{c},x_{2}^{c},\cdots,x_{d-1}^{c}}\sigma_{d,\vec{x}^{c}}\otimes\rho_{x_{1}^{c}}\otimes\rho_{x_{2}^{c}}\otimes\cdots\otimes\rho_{x_{d-1}^{c}},$
(52)

because both $x_{k}$ and $x_{k}^{c}$ sum over the same range, namely
$\\{0,1\\}$. Now it is clear that this term is indeed invariant by a simple
change of dummy variables.

We end this section by noting the peculiarity of the situation: through its
dependence on $S_{\mu\nu}$, the magnetic inversion symmetry $A_{d}$ depends on
the disorder realization of the ensemble, hence the symmetry itself is
disordered. This is exceptional in that the symmetries of most disordered
systems do not depend on disorder realizations. However, the effects of this
disordered symmetry are as real as the conventional cases. In particular, to
study the level statistics we must focus on one block of the $H_{d}$ that is
irreducible under $A_{d}$.

##  5 Sum rules for the Hamiltonian

There are exact sum rules for the Hamiltonian that are valid even without
taking the disorder average. They will account for some salient features of
the level statistics we are going to see in section 6. The sum rules are
consequences of the hypercubic geometry and the fact that the Hamiltonian (3)
has only nearest neighbor hoppings. The sum rules can be evaluated in the
tensor product representation (28) as well. Since this calculation for
$\textrm{Tr}H^{4}$ and $\textrm{Tr}H^{6}$ in tensor product representation is
rather lengthy, we have moved it to appendix A.

###  5.1 $\textrm{Tr}AH^{p}$

In the study of level statistics, we analyze the energy eigenvalues in the
same block under magnetic inversion symmetry $A$. So instead of the total
moments $\textrm{Tr}H^{p}$, what we really should be interested in is
$\frac{1}{2}\textrm{Tr}\left((\mathbb{1}+A)H^{p}\right)$. However in this
section we will see

$\textrm{Tr}A_{d}H_{d}^{p}=0\quad\text{ for}\ p<d.$ (53)

So for low moments we might as well just study $\textrm{Tr}H^{p}$.

#### Geometric picture

The magnetic inversion $A_{d}$ has the physical meaning of particle hopping
from one lattice site to the site sitting on the corresponding longest
diagonal. The $H^{p}$ involves $p$-step hoppings connecting nearest neighbors.
For the trace to be nonzero, we must have at least one hopping configuration
that forms a loop. This means some of the $p$-step nearest-neighbor hops must
reach the longest diagonal to form a loop with the $A_{d}$ hopping. This is
clearly impossible for $p<d$.

#### Tensor product picture

Since $A_{d}$ is the tensor product of $d$ off-diagonal Pauli matrices (see
equation (49)), while each factor contributing the Hamiltonian (28) contains
only one off-diagonal Pauli matrix. It is clear that
$\textrm{Tr}A_{d}H_{d}^{p}=0$ for $p<d$, because for $p<d$ every term in
$H^{p}$ will have at least one diagonal $2\times 2$ matrix in the tensor
product.

###  5.2 $\textrm{Tr}H^{2}$

We wish to prove

$\displaystyle 2^{-d}\textrm{Tr}H_{d}^{2}=d.$ (54)

In fact, we will prove a stronger identity for the diagonal entries of
$H^{2}$:

$\left(H_{d}^{2}\right)_{ii}=d.$ (55)

#### Geometric picture

We note the diagonal elements $(H^{2})_{ii}$ only receive contributions from
2-step loops. But a 2-step loop must be one step through some lattice link
followed by one step back through the same link, and hence the phases cancel.
We can choose to first step to be along any direction, thus in $d$ dimensions
we have $d$ contributions, each being $1$. This gives (55).

#### Tensor product picture

The tensor products in the Hamiltonian (28) involve both diagonal and off-
diagonal two by two matrices. To contribute to $(H_{d}^{2})_{ii}$, terms with
the off-diagonal Pauli matrix must be in the same position in both factors of
$H$, so

$\displaystyle(H_{d}^{2})_{ii}=\left(\sigma_{0}\otimes
H_{d-1}\right)^{2}_{ii}+\left(\sum_{x_{1},x_{2},\cdots,x_{d-1}}\sigma_{d,\vec{x}}\otimes\rho_{x_{1}}\otimes\rho_{x_{2}}\otimes\cdots\otimes\rho_{x_{d-1}}\right)^{2}_{ii},$
(56)

and inside the second term, we have terms

$\displaystyle\left(\sigma_{d,\vec{x}}\otimes\rho_{x_{1}}\otimes\cdots\otimes\rho_{x_{d-1}}\right)\cdot\left(\sigma_{d,\vec{x}^{\prime}}\otimes\rho_{x_{1}^{\prime}}\otimes\cdots\otimes\rho_{x_{d-1}^{\prime}}\right),$
(57)

which are only nonzero if $x_{k}=x_{k}^{\prime}$ for all $1\leq k\leq d-1$. It
is not hard to see then the sum over $x_{1},\ldots,x_{d-1}$ results in a
tensor product of $d$ identity matrices. The same argument applies to
$\left(\sigma_{0}\otimes H_{d-1}\right)$ through the recursive definition of
$H_{d-1}$. We can do this recursively all the way to $H_{1}$ in $d-1$ steps,
and each step creates an identity matrix, so

$(H_{d}^{2})_{ii}=\left(\sigma_{0}\otimes\sigma_{0}\otimes\cdots\otimes
H_{1}\right)^{2}_{ii}+(d-1)=d.$ (58)

###  5.3 $H^{2}$ at $\pi$ flux

We just demonstrated that the diagonal entries of $(H_{d}^{2})_{ii}=d$. We
shall further show that at $\pi$ flux,

$H_{d}^{2}(\phi=\pi)=d\mathbb{1}.$ (59)

Together with the sublattice symmetry described in section 4.1, this implies
at $\pi$ flux $H_{d}$ has exactly half of its eigenvalues being $-\sqrt{d}$
and and the other half being $\sqrt{d}$.

#### Geometric picture

We only need to show

$\left(H_{d}^{2}\right)_{ij}=0\quad\text{if $i\neq j$.}$ (60)

Note $\left(H_{d}^{2}\right)_{ij}$ receives contributions from 2-step lattice
paths that connect lattice sites $i$ and $j$.222Perhaps it is more precise to
say “sites represented by $i$ and $j$”, namely sites whose coordinates are
$\vec{x}(i)$ and $\vec{x}(j)$, whose components are the binary digits of $i-1$
and $j-1$ in reverse order. There are two scenarios for $i\neq j$:

  1. 1.

There is no 2-step path from $i$ to $j$. For such pairs of $ij$ clearly
$(H^{2}_{d})_{ij}=0$.

  2. 2.

Sites $i$ and $j$ can be connected by a 2-step path. If so then sites $i$ and
$j$ must be siting on the diagonal of a face of the hypercube and there are
exactly two paths connecting them, which form the four sides of the face, see
figure 1. If the direction of one of the two paths is reversed, we will have a
Wilson loop of flux $\pi$, and this means the two original paths give
contributions that differ by a factor of $e^{i\pi}=-1$, so their sum vanishes.

Hence (60) is proven.

$i$$j$$\pi$ Figure 1: Two lattice paths that connect sites $i$ and $j$ with a
$\pi$ flux.

#### Tensor product picture

We have seen in the last section that the diagonal entries of $H^{2}$ come
from individual terms squared. Now we need to show the cross terms cancel out
for $\phi=\pi$. One such pair of cross terms is an anticommutator

$\left\\{\sum_{x_{1}^{\prime},\cdots,x_{d-2}^{\prime}}\sigma_{0}\otimes\sigma_{d-1,\vec{x}^{\prime}}\otimes\rho_{x_{1}^{\prime}}\otimes\cdots\otimes\rho_{x_{d-2}^{\prime}},\sum_{x_{1},\cdots,x_{d-1}}\sigma_{d,\vec{x}}\otimes\rho_{x_{1}}\otimes\rho_{x_{2}}\otimes\cdots\otimes\rho_{x_{d-1}}\right\\}$
(61)

where the first factor is part of $\sigma_{0}\otimes H_{d-1}$ and the second
factor is the second term in equation (28). Since $\rho_{0}\rho_{1}=0$, the
product is only nonzero when

$x_{1}^{\prime}=x_{2},\ x_{2}^{\prime}=x_{3},\ldots,x_{d-2}^{\prime}=x_{d-1}.$
(62)

So the sum reduces to333We remind the readers that
$\sigma_{d-1,(x_{2},x_{3},\ldots,x_{d-1},x_{d-1}^{\prime})}$ does not depend
on the last coordinate $x_{d-1}^{\prime}$.

$\displaystyle\sum_{x_{1},x_{2},\cdots,x_{d-1}}\sigma_{d,\vec{x}}\otimes\left\\{\sigma_{d-1,(x_{2},x_{3},\ldots,x_{d-1},x_{d-1}^{\prime})},\rho_{x_{1}}\right\\}\otimes\rho_{x_{2}}\otimes\cdots\otimes\rho_{x_{d-1}}$
(63) $\displaystyle=$
$\displaystyle\sum_{x_{1},x_{2},\cdots,x_{d-1}}\sigma_{d,\vec{x}}\otimes\sigma_{d-1,(x_{2},x_{3},\ldots,x_{d-1},x_{d-1}^{\prime})}\otimes\rho_{x_{2}}\otimes\cdots\otimes\rho_{x_{d-1}},$
(64)

where we have used

$\displaystyle\sigma_{d-1,(x_{2},x_{3},\ldots,x_{d-1},x_{d-1}^{\prime})}\rho_{x_{1}}$
$\displaystyle=\rho_{x_{1}^{c}}\sigma_{d-1,(x_{2},x_{3},\ldots,x_{d-1},x_{d-1}^{\prime})},$
(65) $\displaystyle\rho_{x_{1}}+\rho_{x_{1}^{c}}$
$\displaystyle=\mathbb{1}_{2\times 2}.$ (66)

Note only the first tensor factor in (64) depends on $x_{1}$. Now we see the
sum over $x_{1}$ already gives zero because

$\sigma_{d,(x_{1},x_{2},\ldots,x_{d})}=-\sigma_{d,(x_{1}^{c},x_{2},\ldots,x_{d})}\quad\text{when
}\phi=\pi.$ (67)

The same argument can be applied to all other mixed terms. For $\phi=\pi$ we
thus demonstrated $H_{d}^{2}=d\mathbb{1}$.

###  5.4 $\textrm{Tr}H^{4}$

We wish to prove $\textrm{Tr}H^{4}$ does not depend on disorder realizations
of $S_{\mu\nu}$.

#### Geometric picture

We need to consider all the 4-step loops on the hypercube. If the path is
backtracking then the loop has zero area, so quite trivially they do not
depend on flux realizations. The only other possibility for a 4-step loop is a
one that travels the four sides of a hypercube face, and its contribution to
the trace is its Wilson loop value $e^{i\phi S_{\mu\nu}}$. However, since each
clockwise loop is accompanied by its counterclockwise counterpart, the
contributions must be functions of $\cos\left(\phi
S_{\mu\nu}\right)=\cos\phi$.444The crucial point is that a 4-step loop can at
most loop around one face of the hypercube. For larger loops when several
faces can be looped around, we generically have
$\cos\left[\phi(S_{\mu_{1}\nu_{1}}+S_{\mu_{2}\nu_{2}}+\cdots)\right]$. We see
in both cases the contributions do not depend on the disorder realization of
$S_{\mu\nu}$.

#### Tensor product picture

The fourth moment can also be worked out in the tensor representation of the
Hamiltonian, see appendix A.1. This allows us to obtain the exact result for
the fourth moment which is in agreement with that obtained in Parisi:1994jg .

###  5.5 $\textrm{Tr}H^{6}$

In this section we prove that $\textrm{Tr}H^{6}$ does not depend on disorder
realizations of $S_{\mu\nu}$.

#### Geometric picture

A six-step loop can at most traverse three different dimensions. Let us first
think about $d=3$. As Parisi argued Parisi:1994jg , in three dimensions the
field strength tensor $\phi S_{\mu\nu}$ can be viewed as a vector, pointing
along one of the longest diagonals of the 3-cube. Hence all possible
realizations of the flux are related to each other by a spatial rotation in
the cubic symmetry group, which implies their Hamiltonians all have the same
spectrum independent of $S_{\mu\nu}$. The loops that contribute to
$\textrm{Tr}H^{6}$ can traverse one, two or three different dimensions. Those
that traverse one and two dimensions are independent of $S_{\mu\nu}$ for
reasons discussed in section 5.4. This implies that for $d=3$ in particular,
the sum of all Wilson loops that traverse three different dimensions is also
independent of realizations of $S_{\mu\nu}$. Now let us consider general $d$.
Since every three different dimensions uniquely define a 3-cube, it is evident
that all loops that traverse three different dimensions can be partitioned
into groups by the 3-cubes they reside in. By the argument just laid out, the
sum of each group of such loops is independent of $S_{\mu\nu}$, and hence the
total sum retains the independence. It is important to separate the
contributions of the loops that traverse three different dimensions from the
rest for this argument to work, because a loop that traverses one or two
dimensions can reside in multiple 3-cubes.

#### Tensor product picture

For the calculation using the tensor representation we also have to
distinguish several cases. Although the calculation is straightforward, the
preponderance of indices makes this calculation rather cumbersome, and we have
moved it to A.2. This calculation shows that the disorder independence of
$\textrm{Tr}H^{6}$ arises because we have just enough terms in the expansion
of $\textrm{Tr}H^{6}$ to cancel the sine-dependent terms of the form $\sin\phi
S_{\mu\nu}$. However, the number of sine-dependent terms grows exponentially
while the number of terms available for canceling sine-dependent terms does
not grow as quickly, so for higher moments we cannot expect disorder
independence. As it turns out the same calculation already fails for
$\textrm{Tr}H^{8}$.

##  6 Chaos on the hypercube

For $\phi=0$ the model is integrable, and has a degenerate spectrum (33). The
degeneracies are lifted at nonzero $\phi$, but the eigenvalues will eventually
flow to $\pm\sqrt{d}$ at $\phi=\pi$, as predicted by equation (60) and the
sublattice symmetry. A figure of the spectral flow as a function of $\phi$ is
shown in figure 2 with the quantum number of the magnetic inversion symmetry
equal to $s=1$ in the left figure and $s=-1$ in the right figure. The flow is
similar to the one of the Maldacena-Qi model. At $\phi=0$ the spectrum and
degeneracies are the same as for the Maldacena-Qi model at infinite coupling.
The degeneracies are lifted at nonzero $\phi$, and at $\phi=\pi/2$ the
spectrum splits into two bands, a feature that is not present in the MQ model.
The ground state of the model is separated from the rest of the spectrum by a
gap, and our numerical results suggest that the gap likely remains finite for
$\phi<\pi/2$ in the thermodynamical limit (see the left figure of figure 2).
We expect that the levels in each subsector become chaotic as soon as the
bands emanating from degenerate eigenvalues start overlapping (at about
$\phi=\pi/4$) which will be studied in more detail below.

Figure 2: Spectral flow of the hypercubic Model as a function of the flux
$\phi$. for each of its two symmetry classes. To make individual curves
visible, we show the results for $d=8$, but the features for larger values of
$d$ are similar.

The apparent crossings of the spectral flow lines are actually avoided
crossings even though some are extremely close.

###  6.1 Average spectral density

It was already realized by Parisi that the spectral density of the large $d$
limit of the hypercube model is given by the spectral density of the
Q-harmonic oscillator. The argument is the same as in the case of the SYK
model Erdos:2014a ; Cotler2016 ; Garcia-Garcia:2017pzl ; Berkooz:2018qkz ;
Berkooz:2018jqr , and can be summarized as follows (see appendix B for more
details). Traces of the lattice Hamiltonian,

$\displaystyle\textrm{Tr}H^{2p},$ (68)

can be written as a sum of Wilson loops on the lattice. As is explained in
appendix B paths can be represented as chord diagrams, with a nonzero Wilson
loop represented by a crossing. Each crossing gives rise to a factor
$q=\cos\phi$.555Note this $q$ is not the $q$ often used in the context of SYK
model where it denotes the interaction order of Majorana fermions. For large
$d$ the leading contributions is from Wilson loops covering the maximum number
dimensions. After ensemble averaging we thus obtain the $2p$-th moment:

$\displaystyle
M_{2p}^{\text{HC}}:=\frac{2^{-d}\langle\textrm{Tr}H^{2p}\rangle}{\left(2^{-d}\langle\textrm{Tr}H^{2}\rangle\right)^{p}}=\sum_{k=0}^{d(d-1)/2}a_{k}q^{k}+O(1/d),$
(69)

where $a_{k}$ is the number of chord diagrams with $k$ crossings. We have
defined $M_{2p}^{\text{HC}}$ as a reduced moment since we used
$\textrm{Tr}H^{2}$ in the denominator, but we will call $M_{2p}^{\text{HC}}$
“moment” when the context is free of confusion. In appendix B we lay out the
arguments and derivations that lead to equation (69) in more details, and
discuss the subleading corrections.

The moments given in equation (69) are the moments of the density function of
the Q-Hermite polynomials:

$\displaystyle\rho^{\rm
QH}(x)=\left(1-\frac{x^{2}}{e_{0}^{2}}\right)^{1/2}\frac{2}{\pi
e_{0}}\prod_{k=1}^{\infty}\frac{1-Q^{2k}}{1-Q^{2k-1}}\prod_{k=1}^{\infty}\left(1-4\frac{x^{2}}{e_{0}^{2}}\frac{Q^{k}}{(1+Q^{k})^{2}}\right),$
(70)

with $e_{0}=2/\sqrt{1-Q}$ and $Q=q=\cos\phi$. However, to include some of the
finite-$d$ corrections we set $Q=\eta$, which is a renormalized version of
$q$, obtained by matching the fourth moment of $\rho^{\text{QH}}(x)$ and the
fourth moment of the hypercube model exactly:

$\displaystyle Q=\eta:=M_{4}^{\rm HC}-2=\cos\phi-\frac{\cos\phi+1}{d}.$ (71)

It is clear $\eta\to q=\cos\phi$ in the large $d$ limit. In figure 3 we show
the average spectral densities for three different values of $\phi$ and
compare the result with the Q-Hermite spectral density with $Q=\eta$.
Renormalizing $q$ to $\eta$ improves the accuracy for finite $d$, but this is
still not exact: the deviation will start to appear for the sixth and higher
moments. We cite Marinari:1995jwr  here for the exact results up to the eighth
moment:

$\displaystyle M_{4}^{\text{HC}}$
$\displaystyle=\frac{d-1}{d}(2+q)+\frac{1}{d},$ (72) $\displaystyle
M_{6}^{\text{HC}}$
$\displaystyle=\frac{(d-1)(d-2)}{d^{2}}(5+6q+3q^{2}+q^{3})+\frac{d-1}{d^{2}}(9+6q)+\frac{1}{d^{2}},$
(73) $\displaystyle M_{8}^{\text{HC}}$
$\displaystyle=\frac{(d-1)(d-2)(d-3)}{d^{3}}(14+28q+28q^{2}+20q^{3}+10q^{4}+4q^{5}+q^{6})$
(74)
$\displaystyle\quad+\frac{(d-2)(d-1)}{d^{2}}(56+86q+52q^{2}+16q^{3})+\frac{d-1}{d^{3}}(33+28q+2q^{2})+\frac{1}{d^{3}}.$

Figure 3: The spectral density of the hyper cubic model (black curve) compared
to the eighth order Q-Hermite spectral density defined in (75) for three
different flux values as given in the legend of the figures. In the left
figure $\phi=0.57\pi/2$, the curve resulted from a fitted $a_{8}$ (red curve)
differ from the one from the $a_{8}$ that is calculated by equation (80) (blue
curve) by less than the line width of the curve’s plot, so we do not see the
red curve at all. The red curve and blue curve also agree very well in the
middle figure where $\phi=0.95\pi/2$. The deviations are larger for
$\phi=\frac{3}{4}\pi$ in the right figure.

In terms of $\rho^{\rm QH}_{\eta}(x)$ the spectral density (before ensemble
averaging) can be expanded as

$\displaystyle\rho^{\rm HC}(x)=\rho^{\rm
QH}_{\eta}(x)(1+a_{6}H^{\eta}_{6}(x)+a_{8}H^{\eta}_{8}(x)+\cdots),$ (75)

where the coefficients $a_{8},\ldots$ are random variables (note that $a_{6}$
is determined by the sixth moment and does not depend on the disorder
realization), and $H^{\eta}_{n}$ are the Q-Hermite polynomials defined by the
recursion relation Viennot-1987

$\displaystyle
H_{n+1}^{\eta}(x)=xH_{n}^{\eta}(x)-\sum_{k=0}^{n-1}\eta^{k}H_{n-1}^{\eta}(x)$
(76)

with the initial conditions

$\displaystyle H_{0}^{\eta}(x)=1\qquad{\rm and}\qquad H_{1}^{\eta}(x)=x.$ (77)

The Q-Hermite polynomials satisfy the orthogonality relation

$\displaystyle\int_{-\frac{2}{\sqrt{1-\eta}}}^{\frac{2}{\sqrt{1-\eta}}}dx\rho^{\text{QH}}_{\eta}(x)H_{n}^{\eta}\left(x\right)H_{m}^{\eta}\left(x\right)=\delta_{nm}n_{\eta}!,$
(78)

where $n_{\eta}!$ is the Q-factorial defined as

$\displaystyle n_{\eta}!=\prod_{k=1}^{n-1}\
\left(\sum_{s=0}^{k}\eta^{s}\right).$ (79)

Note that for the choice of $\eta$ in (71) the coefficients of
$H^{\eta}_{2}(x)$ and $H^{\eta}_{4}(x)$ vanish since $\rho^{\text{QH}}_{\eta}$
already gives the exact results for $M^{\text{HC}}_{2}$ and
$M^{\text{HC}}_{4}$. We stress that they vanish not just after averaging but
also realization by realization, this is because in section 5 we have proven
$\textrm{Tr}H^{2}$ and $\textrm{Tr}H^{4}$ are independent of disorder
realizations. The coefficients $a_{6}$ and $a_{8}$ after ensemble averaging
are given by (in the normalization where $M_{2}=1$)

$\displaystyle a_{6}$ $\displaystyle=$
$\displaystyle\frac{1}{6_{\eta}!}(M_{6}^{\rm HC}-M_{6}^{{\rm QH},\eta}),$
$\displaystyle\langle a_{8}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{8_{\eta}!}(M_{8}^{\rm HC}-M_{8}^{{\rm
QH},\eta})-\frac{(7+6\eta+5\eta^{2}+4\eta^{3}+3\eta^{4}+2\eta^{5}+\eta^{6})a_{6}6_{\eta}!}{8_{\eta}!},$
(80)

where we note again $\langle a_{6}\rangle=a_{6}$ because $\textrm{Tr}H^{6}$ is
independent of disorder realizations, which is not true for $a_{8}$. This is
not a good expansion for negative $\eta$ when $n_{\eta}!$ becomes small, see
table 1. For example, the large value of $a_{8}$ for $\phi=3\pi/4$ is due to
the smallness of $8_{\eta}!\approx 0.01$. The expansion diverges for
$\phi\to\pi$. The reason is that

$\displaystyle(2p)_{(-1+x)}!=p!x^{p}+O(x^{p+1}),$ (81)

while

$\displaystyle M_{8}^{\rm HC}-M_{8}^{{\rm
QH},q-(q+1)/d}\sim-\frac{(q+1)^{2}}{3d},$ (82)

so that $a_{8}$ diverges as $1/(q+1)^{2}$ for $q\to-1$. This explains why in
the left two figures of figure 3 the fitted values of $a_{8}$ are close to the
calculated values of $a_{8}$, whereas the in the right figure the agreement is
not as good. For $a_{6}$ we are in a better position. Then

$\displaystyle M_{6}^{\rm HC}-M_{6}^{{\rm
QH},q-(q+1)/d}\sim\frac{(q+1)^{3}}{d^{2}}.$ (83)

so that $a_{6}\sim 1/d^{2}$. This also explains that $a_{6}<<1$ for $d=14$,
see table 1. For a given realization, the expansion coefficient $a_{8}$ is
also given by equation (80) but with $M_{8}^{\text{HC}}$ replaced by the
eighth moment of that realization.

###  6.2 Spectral correlations

Figure 4: The number variance $\Sigma^{2}(n)$ versus $n$ for $\phi=0.57\pi/2$,
$\phi=0.95\pi/2$ and $\phi=\frac{3}{4}\pi$. The spectra have been unfolded
using the ensemble average of the spectral density. The right figure shows the
number variance for larger values of $n$.

In the SYK model spectral correlations show agreement with random matrix
theory for a distance of about $2^{N/2}/N$ level spacings if the fluctuations
from one realization to the next one are eliminated. If we include those
fluctuations, the range of agreement is reduced to $O(N^{2})$ which can be
easily understood by analyzing the effect of overall scale fluctuations due to
the fact that the number of independent random variables is only of order
$N^{4}$ Flores_2001 ; Altland:2017eao ; Garcia-Garcia:2018ruf ;
Gharibyan:2018jrp ; Jia:2019orl  while the number of eigenvalues is
$2^{N/2}/2$. In the hypercubic model, the first six moments are independent of
the realization, and fluctuations of the overall scale and low-order moments
are mostly absent. The sixth order Q-Hermite result already gives a very
accurate description of the average spectral density for values of
$\frac{\pi}{4}<\phi<\pi/2$. Indeed for $\phi=0.57\pi/2$, there is very little
difference in the statistical spectral observables between local unfolding,
where the spectral density of each realization is fitted to a smooth curve,
and unfolding with the ensemble-averaged spectral density. In the left column
of figure 4 we show the number variance $\Sigma^{2}(n)$ versus the number of
levels in an interval containing $n$ levels on average for $n$ up to 50, and
in the right column (black curves) up to 1000. In figure 5 we show the same
quantities but with local unfolding. We compare these results to the
analytical expression for the Gaussian Unitary Ensemble (red curve).
Deviations from the universal random matrix curve start at $n\approx d$. This
is in agreement with the observation that the hypercubic Hamiltonian is
determined by $O(d^{2})$ random variables so that the relative fluctuations in
$a_{8}$ and higher order expansion coefficients are of order $1/d$.

$\cos\phi$ | $\eta$ | $a_{6}$ | $\langle a_{8}\rangle$ | $\langle(\delta a_{8})^{2}\rangle^{1/2}$ | $\langle a_{8}\rangle^{\text{fit}}$
---|---|---|---|---|---
0.6252 | 0.509 | -0.0010 | $-4.16\times 10^{-4}$ | $2.66\times 10^{-4}$ | $-3.72\times 10^{-4}$
1/13 | 0 | -0.0059 | -0.00237 | 0.0273 | $-3.68\times 10^{-4}$
$-1/\sqrt{2}$ | -0.728 | -0.0038 | -2.003 | 0.086 | -1.80
Table 1: Collective spectral fluctuations as measured by the coefficient
$a_{8}$ in the expansion of the spectral density in Q-Hermite polynomials. The
fitted value of $a_{8}$ is within the range of these fluctuations.

The fluctuations of the number of levels in an interval containing $n$ levels
on average is thus $\delta n/n\sim O(1/d)$ resulting in a correction to the
number variance that behaves as $n^{2}/d^{2}$. The results for
$\phi=0.57\pi/2$ are significantly closer to the random matrix result than
those for the other values of $\phi$. For the first ($\phi=0.57\pi/2$) and
second row ($\phi=\arccos(1/13)$) of figure 4 we used the ensemble average of
the eighth order Q-Hermite result to unfold the spectral density, while for
the third row ($\phi=\frac{3}{4}\pi$) a third order polynomial fit to the
ensemble average of the spectral density was used to unfold the bulk of the
spectrum.

Figure 5: The number variance, $\Sigma^{2}(n)$ versus the average number of
eigenvalues in the interval for $d=14$ and values of the flux as indicated in
the legends of the figures. The results have been obtained by unfolding the
spectral for each realization separately (local unfolding). For
$\phi=0.57\pi/2$, the curves are indistinguishable from the results for
ensemble unfolding. The left figures give the same curves as the right figures
but for a smaller range of $n$.

Figure 6: The ensemble-unfolded spectral form factor for $d=14$ at
$\phi=0.57\pi/2$, $\phi=0.95\pi/2$ and $\phi=3\pi/4$ (black curves). The
results are compared with the result for the GUE (red curves). In the left
figure, the careful observer can see a tiny peak at $\tau$ close to zero which
is responsible for the large deviation of the number variance from the
universal GUE result. This peak is magnified in the right figure (black
curve), where we also show the result for local unfolding (blue curve).
Ensemble unfolding and local unfolding give almost indistinguishable results
for $\phi=0.57\pi/2$. If the results of local unfolding were plotted in the
left figures, the differences with the ensemble unfolding results would not be
visible for any of the three figures.

The difference between the results for ensemble unfolding and local unfolding
is due to the fluctuations of $a_{8}=\langle a_{8}\rangle+\delta a_{8}$. Table
1 contains the results for the simulation parameters of the above figures. We
conclude that for $\phi=0.57\pi/2$ the collective fluctuations only contribute
negligible to the spectral fluctuations, while they are important for
$\phi=\arccos(1/13)$ and $\phi=3\pi/4$.

The deviations from the universal RMT result are barely visible in the
spectral form factor (see the left column of figure 6), where the results for
the hypercube model (black curve) agree very well with the GUE result (red
curve) except for a very narrow peak for $t$ close to zero. To reduce finite
size effects, the spectral form factor is calculated using a Gaussian window
of width 2000 for $\phi=0.57\pi/2$ and $\phi=\arccos(1/13)$; for
$\phi=3\pi/4$, where the range of the spectrum where the spectrum can be
reliably unfolded is smaller, the width is taken to be $500$. For
$\phi=0.57\pi/2$ local unfolding and ensemble unfolding give almost identical
results (see upper right figure of figure 6), while for the other values of
$\phi$ in this figure, there are significant reductions of the small time
peaks for local unfolding (blue curves). This suggests the moments that are
responsible for the early-time peak are much beyond the eighth order, and more
so for $\phi=0.57\pi/2$ than larger values of $\phi$. Indeed, as we have shown
in section 5, there is no fluctuation up to the sixth moment, so that the
first moment that can fluctuate is the eighth moment. In this light it is
perhaps not too surprising that the eighth-order local unfolding does not
eliminate much fluctuation. It is instructive to contrast this phenomenon in
the HC model to its counterpart in the SYK model Jia:2019orl , where the
eighth-order local unfolding is quite adequate to remove the early-time peak
that is present in the ensemble-unfolded spectral form factor. The early-time
peak is responsible for the deviation from the random matrix result. This can
be shown explicitly by calculating the number variance directly from the
spectral form factor with and without this peak using the relation delon-1991

$\displaystyle\Sigma^{2}(n)=\frac{n^{2}}{2\pi^{2}}\int_{-\infty}^{\infty}dtK(t)\left(\frac{\sin(nt/2)}{nt/2}\right)^{2}.$
(84)

Note the derivation of this relation assumes translational invariance of the
spectral correlations which is not the case close to the center of the
spectrum for a chirally symmetric spectrum.

Since we deal with a bipartite lattice the Hamiltonian has a chiral symmetry,
and the eigenvalues correlations are in the universality class of chiral
Random Matrix Theory Verbaarschot:1994qf , specifically the chiral Gaussian
Unitary Ensemble (chGUE) since the system does not have any anti-unitary
symmetry. The chGUE ensemble is characterized by an oscillatory structure in
the spectral density near zero on the scale of the average level spacing, and
we call the spectral density in this regime the microscopic spectral density.
The microscopic spectral density is defined by Shuryak:1992pi

$\displaystyle\rho_{s}(E)=\frac{1}{\Sigma N}\rho\left(\frac{E}{\Sigma
N}\right),$ (85)

where

$\displaystyle\Sigma=\lim_{\lambda\to
0}\lim_{N\to\infty}\frac{\pi\rho(\lambda)}{N}$ (86)

and $N$ is a parameter that counts the total number of eigenvalues such as the
size of the random matrix. For an overview of chiral Random Matrix Theory and
its applications to lattice QCD we refer to Verbaarschot:2000dy . In the case
of hypercube model $N=2^{d}$ and
$\rho(\lambda)=\langle\rho^{HC}(\lambda)\rangle$. In figure 7 we show the
microscopic spectral density for an ensemble of 10,000 Hamiltonians for $d=12$
and $\phi=0.57\pi/2$ (black dots). The result is compared with the analytical
result for the chGUE microscopic spectral density (red curve)
Verbaarschot:1993pm :

$\displaystyle\rho_{s}(E)=\frac{E}{2}(J_{0}^{2}(E)+J_{1}^{2}(E)),$ (87)

where $J_{n}(E)$ are the Bessel functions. We remark that there is no fitting
and the agreement is excellent.

Figure 7: The microscopic spectral densities. Black dots: numerically
calculated from 10,000 realizations of the $d=12,\phi=0.57\pi/2$ Hamiltonian.
Red curve: analytically predicted by the chGUE random matrix theory.

The chiral symmetry also affects the number variance, but the effects are
negligible unless the intervals for which the number variance is calculated
are chosen symmetrically about zero. The correlations due to the pairing
$\pm\lambda_{k}$ are also visible in the short time behavior of the form
factor. Instead of $K_{c}(t)\sim t^{2}$ for the GUE we have $K_{c}(t)\sim
t^{4}$ for the chGUE, both for $t\to 0$ at finite size of the matrices.
However, the peak near zero in the numerical results obscures this effect. The
number variance of the chGUE is reduced by a factor 2 (in the domain where
$\Sigma^{2}(n)\sim\log n$) for intervals that are symmetric about zero
Toublan:2000dn . However, because we calculate the number variance by spectral
averaging over the spectrum, this has only a small effect except when $n$
becomes large. In fact the kinks in the number variance for $n>400$ are due to
this effect.

##  7 Thermofield double state

In this section we construct the ThermoField Double (TFD) state corresponding
the ground state of the hypercubic model. Whether or not this is the case is a
basis-dependent statement, and we have to identify an appropriate basis.
Inspired by the Maldacena-Qi model we use the sum of a left SYK model and a
right SYK model, and in this case we illustrate our construction by choosing a
two-body Hamiltonian. We remark that in the MQ model, “left” and “right” refer
to the two sides of a worm hole, and quantum mechanically this translates to
the fact the elementary fermion operators factorize into tensor products in a
product Hilbert space. In this paper we do not dwell on the space-time
interpretations of the HC model, so we use the terms to simply refer to the
tensor product structure. General arguments to construct a TFD state were
given is cottrell:2018ash , and applications of the TFD state can be found in
Maldacena:2001kr ; delCampo:2017ftn .

The first observation is that the coupling of the Maldacena-Qi model is
equivalent to the Parisi Hamiltonian at zero flux, which can be expressed in
terms of the gamma matrices defined in equation (31). We thus have

$\displaystyle
H(\phi=0)=i\sum_{k=1}^{d}\gamma_{k}^{L}\gamma_{k}^{R}=UH_{MQ}U^{-1}$ (88)

with

$\displaystyle
H_{MQ}=i\sum_{k=1}^{d}\tilde{\gamma}_{k}^{L}\tilde{\gamma}_{k}^{R},$ (89)

where the gamma matrices $\tilde{\gamma}_{k}^{L(R)}$ are in a representation
that was used in Garcia-Garcia:2019poj  to prove that the ground state of the
Maldacena-Qi model is a TFD state. Specifically,

$\displaystyle\begin{split}\tilde{\gamma}_{k}^{L}=&\tilde{\gamma}_{k}\otimes
1,\quad k=1,2,\ldots,d/2,\\\
\tilde{\gamma}_{k}^{R}=&\tilde{\gamma}_{c}\otimes\tilde{\gamma}_{k},\quad
k=1,2,\ldots,d/2,\end{split}$ (90)

where $\tilde{\gamma}_{k}$ are Dirac matrices in $d/2$ dimensions and
$\tilde{\gamma}_{c}$ is the corresponding chirality Dirac matrix. For this
construction to work we need $d/2$ to be even, namely $d$ is a multiple of
$4$. The $\tilde{\gamma}^{L}$ and $\tilde{\gamma}^{R}$ matrices can be
obtained by a permutation of the $\gamma^{L}$ and $\gamma^{R}$ matrices in
equation (31) as follows:

$\begin{split}\tilde{\gamma}^{L}_{2k-1}&=\gamma^{L}_{k},\
\tilde{\gamma}^{L}_{2k}=\gamma^{R}_{k},\\\
\tilde{\gamma}^{R}_{2k-1}&=\gamma^{L}_{\frac{d}{4}+k},\
\tilde{\gamma}^{R}_{2k}=\gamma^{R}_{\frac{d}{4}+k},\end{split}$ (91)

for $k=1,2,\ldots,d/4$. Then we can check the $\tilde{\gamma}_{k}$ matrices in
equation (90) take the form:

$\begin{split}\tilde{\gamma}_{2k-1}&=\overbrace{\sigma_{1}\otimes\cdots\otimes\sigma_{1}}^{k}\otimes\sigma_{3}\otimes\overbrace{\sigma_{0}\otimes\cdots\otimes\sigma_{0}}^{\frac{d}{2}-k-1},\\\
\tilde{\gamma}_{2k}&=\overbrace{\sigma_{1}\otimes\cdots\otimes\sigma_{1}}^{k}\otimes\sigma_{2}\otimes\overbrace{\sigma_{0}\otimes\cdots\otimes\sigma_{0}}^{\frac{d}{2}-k-1},\end{split}$
(92)

for $k=1,2,\ldots,d/4$, and

$\displaystyle\tilde{\gamma}^{c}=\overbrace{\sigma_{1}\otimes\cdots\otimes\sigma_{1}}^{\frac{d}{2}}.$
(93)

Since both $\\{\gamma^{L}_{k},\gamma^{R}_{k}\\}$ and
$\\{\tilde{\gamma}^{L}_{k},\tilde{\gamma}^{R}_{k}\\}$ are Hermitian
representations of the Clifford algebra in even dimensions, the similarity
transformation $U$ in equation (88) that relates the two is unitary. In the
Maldacena-Qi model, the basis of the TFD state is constructed from the
Hamiltonian

$\displaystyle
H_{SYK}^{R}+H_{SYK}^{L}=\sum_{\alpha}J_{\alpha}(\tilde{\Gamma}_{\alpha}^{R}+\tilde{\Gamma}_{\alpha}^{L}),$
(94)

where $\tilde{\Gamma}_{\alpha}^{L(R)}$ is a product of $q$ different
$\tilde{\gamma}^{L(R)}$ matrices, $\alpha$ is the set of $q$ indices of these
gamma matrices, and $J_{\alpha}$ is the Gaussian-random coupling.666Note again
that here $q$ is an integer in the SYK model, independent of the HC model’s
flux parameter. It is important that the left and right Hamiltonian share the
same coupling $J_{\alpha}$. Because of the tensor structure of the Hamiltonian
it is clear that the eigenstates of this Hamiltonian are given by

$\displaystyle|m\rangle\otimes|n\rangle$ (95)

with eigenvalues $E_{m}+E_{n}$. Here, $|m\rangle$ are the eigenstates of
$H^{L}$ projected onto the the left space. In this basis, the thermofield
double state at inverse temperature $\beta$ is given by

$\displaystyle|{\rm TDF}\rangle=2^{-d/4}\sum_{m}e^{-\beta
E_{m}}|m\rangle|C^{R}e^{\frac{\pi}{4}i\tilde{\gamma}_{c}}Km\rangle$ (96)

with $C^{R}$ the charge conjugation matrix,

$\displaystyle{C^{R}}^{\dagger}\tilde{\gamma}_{k}^{R}C^{R}=\tilde{{\gamma}_{k}^{R}}^{*},$
(97)

and $K$ the complex conjugation operator. In a convention where gamma matrices
$\tilde{\gamma}_{2k}$ are purely imaginary while the $\tilde{\gamma}_{2k-1}$
are purely real like in equation (92), we have that

$\displaystyle C^{L(R)}=\prod_{k=1}^{d/4}{\tilde{\gamma}_{2k-1}^{L(R)}}.$ (98)

The argument to show that the ground state of the Hamiltonian $H_{MQ}$ is
given by the TFD state at $\beta=0$ does not depend on the details of the
Hamiltonian (94) that determines the basis states Garcia-Garcia:2019poj , for
example it does not matter if we use a 2-body, 4-body or 6-body SYK model
Hamiltonians. This follows from the expectation value

$\displaystyle 2^{-d/2}\sum_{mn}\langle m|\langle
C^{R}e^{\frac{\pi}{4}i\tilde{\gamma}_{c}}Km|i\sum_{k}\tilde{\gamma}_{k}^{L}\otimes\tilde{\gamma}_{k}^{R}|n\rangle|C^{R}e^{\frac{\pi}{4}i\tilde{\gamma}_{c}}Kn\rangle$
(99) $\displaystyle=$ $\displaystyle 2^{-d/2}\sum_{mn}\langle
m|{\tilde{\gamma}_{k}^{L}}\tilde{\gamma}_{c}|n\rangle\langle
Km|e^{-\frac{\pi}{4}i\tilde{\gamma}_{c}}{C^{R}}^{\dagger}i{\tilde{\gamma}_{k}^{k}}C^{R}e^{\frac{\pi}{4}i\tilde{\gamma}_{c}}|Kn\rangle$
$\displaystyle=$ $\displaystyle 2^{-d/2}\sum_{mn}\langle
m|{\tilde{\gamma}_{k}}\tilde{\gamma}_{c}|n\rangle\langle
Km|\tilde{\gamma}_{c}{{\tilde{\gamma}_{k}}}^{*}|Kn\rangle$ $\displaystyle=$
$\displaystyle 2^{-d/2}\sum_{mn}\langle
m|{\tilde{\gamma}_{k}}\tilde{\gamma}_{c}|n\rangle\langle
n|{\tilde{\gamma}_{k}}\tilde{\gamma}_{c}|m\rangle,$

where in going from the second line to the third line, we have used the fact
that

$\tilde{\gamma}_{c}^{*}=\tilde{\gamma}_{c},\quad
ie^{-i\frac{\pi}{2}\tilde{\gamma}_{c}}=\tilde{\gamma}_{c}.$ (100)

Now we can use completeness to do the sum over $n$, and after employing that
the gamma matrices square to 1, the sum over $k$ yields a factor $d$ resulting
in

$\displaystyle\langle\text{TFD}|i\sum_{k}\tilde{\gamma}_{k}^{L}\otimes\tilde{\gamma}_{k}^{R}|\text{TFD}\rangle=-d.$
(101)

Since $-d$ is the ground state energy and the ground state is nondegenerate,
the TFD state must be the ground state.

To illustrate the above argument, we choose the two-body SYK Hamiltonian

$\displaystyle
H_{SYK}=\sum_{k<l}J_{kl}\left(i\gamma_{k}^{L}\gamma_{l}^{L}+i\gamma_{k}^{R}\gamma_{l}^{R}\right),$
(102)

to determine the basis states entering the the TDF state and consider the
overlap with the ground state of

$\displaystyle H(\phi=0)=i\sum_{k}\gamma_{k}^{L}\gamma_{k}^{R}.$ (103)

The gamma matrices in both Hamiltonians are in the representation (31). Since
the overlap between states is invariant under a unitary transformation, we can
do the unitary transformation $U$ in equation (88) to transform the
Hamiltonians (102) and (103) into the Hamiltonians (94) and the coupling
matrix in the right-hand side of (89), respectively. Using the above argument,
the ground state of (103) is given by

$\displaystyle
U^{-1}|\text{TFD}\rangle=U^{-1}2^{-d/4}\sum_{m}|m\rangle|C^{R}e^{\frac{\pi}{4}i\tilde{\gamma}_{c}}Km\rangle.$
(104)

Since for even $d/2$ the anti-commutator $\\{C^{R}K,H^{R}_{SYK}\\}=0$, if
$|m\rangle$ is an eigenstate of $H^{R}_{SYK}$ with eigenvalue $E_{m}$, then
$C^{R}K|m\rangle$ is an eigenstate of $H^{R}_{SYK}$ with eigenvalue $-E_{m}$.
The ground state of (103) is thus a linear combination of the zero energy
states of (102). In figure 8 we show the magnitude of the overlap of the
ground state with the $|m\rangle|C^{R}e^{\frac{\pi}{4}i\gamma_{c}}Km\rangle$
(denoted by $|m\rangle|-m\rangle$ in the figure) for $d=12$. The total
strength in this subspace decreases rapidly with increasing magnetic flux, but
the temperature of the TFD state remains infinite.

Figure 8: The overlap of the ground state of the HC model with the components
of the TFD state (left) for several values of the magnetic flux. In the right
figure we show the total overlap between the ground state of the HC model and
the TFD state. Because the zero energy states (102) are degenerate we have a
small symmetry breaking term to this Hamiltonian by
$i\gamma_{k}^{L}\gamma_{l}^{L}\to i\gamma_{k}^{L}\gamma_{l}^{L}(1+\epsilon)$.
For the figures above we used $\epsilon=10^{-4}$.

There are other possibilities to choose a basis for a TFD state. For example
at zero flux, the Hamiltonian may be written as

$\displaystyle H_{d}(\phi=0)=H_{d/2}(\phi=0)\otimes 1+1\otimes
H_{d/2}(\phi=0),$ (105)

and a TFD state can be constructed out of the eigenstates of $H_{d/2}$. For
$\phi\neq 0$ the Hamiltonian

$\displaystyle H_{d/2}(\phi)\otimes 1+1\otimes H_{d/2}(\phi)$ (106)

has its nonzero matrix elements in the same position as the ones of
$H_{d}(\phi)$, and also at $\phi\neq 0$ the eigenstates of $H_{d/2}$ could be
used to construct a TFD state. We have explored these and other related
possibilities, but they did not give a better description of the ground state
of the hypercubic Hamiltonian.

##  8 Conclusions

We have studied the spectral density and the spectral correlations of Parisi’s
hypercubic model. This model is described by the Laplacian on a hypercube with
two lattice points in each dimension and U(1) gauge fields on the links such
that the magnitude of the magnetic flux through each of its faces is constant,
but its orientation is chosen to be random. We have confirmed that the
spectral density of this model is given by the density function of the
Q-Hermite polynomials. This has the important implication that the spectral
density above the ground state $E_{0}$ behaves as $\sinh\sqrt{(E-E_{0})}$.
However, contrary to the SYK model, the ground state is separated from the
rest of the spectrum by a gap. In this respect, the hypercubic model resembles
the Maldacena-Qi model, and we expect it to have a similar phase diagram with
a first order phase transition as a function of the temperature. We hope to
address this point in a future publication. Remarkably, at zero flux the
Hamiltonian of this model coincides with the coupling term of the Maldacena-Qi
model. We have constructed a basis such that in the zero-flux case the ground
state is given by a thermofield double state. Contrary to the Maldacena-Qi
model, at nonzero flux the overlap with the TFD state rapidly decreases. Since
the hypercubic Hamiltonian at nonzero flux is not the sum of a left and a
right Hamiltonian, this did not come as a surprise.

The initial analysis of the spectral correlations of this model led to the
observation that they are described by the superposition of two Gaussian
Unitary Ensembles. This resulted in the discovery of a discrete symmetry which
we later identified as a magnetic inversion symmetry which are analogous to
magnetic translation symmetries studied in the literature. Since this operator
is related to space inversion (which is the same as a translation mod 2 on a
hypercubic lattice), it squares to unity and its eigenvalues are $\pm 1$. We
have analyzed the correlations of the eigenvalues of the hypercubic
Hamiltonian for fixed quantum number of this symmetry and found that they are
correlated according to the GUE. Since this model is determined by $d^{2}$
random numbers, the fluctuations of the number of eigenvalues in an interval
containing $n$ eigenvalues on average behave as $\delta n/n\sim 1/d$, and
hence the number variance for large $n$ behaves as $\Sigma^{2}(n)\sim
n^{2}/d^{2}$ resulting in a “Thouless energy scale” of order $d$. This is in
qualitative agreement with our results. In the spectral form factor, this
deviation is visible as a peak close to zero with area $\sim 1/d^{2}$, which
is only apparent in plots of the connected form factor (which we always plot).

Because of the sublattice symmetry, the Hamiltonian has a chiral symmetry with
eigenvalues occurring in pairs $\pm\lambda_{k}$ so that the eigenvalues are
correlated according the the chiral Gaussian Unitary Ensemble (chGUE). Indeed
we have shown that the microscopic spectral density shows the universal
oscillations characteristic for this ensemble. If the number variance is
calculated for interval that is symmetric about zero, the chiral symmetry
reduces the variance by a factor two. Since we calculate the number variance
by averaging over the spectrum, this effect only affects large values of $n$
where the number variance is dominated by the $n^{2}$ correction.

The traces of powers of the hypercubic Hamiltonian are given by the Wilson
loops of closed paths on the hypercube. We have extended (in appendix B)
Parisi’s work on a one-to-one mapping between these paths and the chord
diagrams that occur in the the calculation of the moments of the SYK
Hamiltonian. This result explains that in both cases the spectral density of
the Hamiltonian is given by the density function of the Q-Hermitian
polynomials. This suggests that the low-energy effective partition function of
the hypercubic Hamiltonian can also be expressed in terms of a Schwarzian
action. We hope to address this point in a future publication.

Our work confirms the strength of random matrix universality. Although the
model is very different from a random matrix theory, and describes a many-body
theory with a sparse Hamiltonian, the level correlations are still very well
described by the corresponding Random Matrix Theory. This further supports the
paradigm, going back to the first applications of random matrix theory to the
nuclear many-body problems, that generically spectra of many-body systems are
chaotic.

###### Acknowledgements.

We acknowledge partial support from U.S. DOE Grant No. DE-FAG-88FR40388. We
acknowledge useful comments by Antonio García-García and Dario Rosa, and their
critical reading of the manuscript. We also acknowledge Peter van
Nieuwenhuizen for a useful discussion.

##  Appendix A Disorder independence of the fourth and sixth moments in the
tensor product representation

###  A.1 $\textrm{Tr}H^{4}$

In this section we calculate $\textrm{Tr}H^{4}$ in the tensor representation
of the Hamiltonian. We obtain an explicit expression for the fourth moment. In
agreement with the geometric picture in the main text, it only depends on the
magnitude of the magnetic flux through the faces of the hypercube and is
independent of its random orientations.

To facilitate the discussion, we define (Here, $\sigma_{0}$ is the $2\times 2$
identity matrix, and we refer to equation (20) for the definition $\rho_{k}$.)

$h_{d,\mu}:=\overbrace{\sigma_{0}\otimes\sigma_{0}\cdots\otimes\sigma_{0}}^{d-\mu}\otimes\sum_{x_{1},\ldots,x_{\mu-1}}\sigma_{\mu,\vec{x}}\otimes\rho_{x_{1}}\cdots\otimes\rho_{x_{\mu-1}},$
(107)

so that the Hamiltonian of the hypercubic model is given by (see equation
(28))

$H_{d}=\sum_{\mu=1}^{d}h_{d,\mu}$ (108)

The fourth moment can be expressed as

$\textrm{Tr}H_{d}^{4}=\sum_{\mu\nu\kappa\omega}\textrm{Tr}\left(h_{d,\mu}h_{d,\nu}h_{d,\kappa}h_{d,\omega}\right).$
(109)

Since each $h_{d,\mu}$ has only one off-diagonal $2\times 2$ matrix in the
tensor product, and its position is labeled by $\mu$, the only nonzero traces
are of the forms $\textrm{Tr}h_{d,\mu}^{4}$,
$\textrm{Tr}\left(h_{d,\mu}h_{d,\mu}h_{d,\nu}h_{d,\nu}\right)$ with $\mu>\nu$,
$\textrm{Tr}\left(h_{d,\mu}h_{d,\nu}h_{d,\mu}h_{d,\nu}\right)$ with $\mu>\nu$,
and their cyclic permutations. It is clear that

$\displaystyle\textrm{Tr}\;h_{d,\mu}^{4}=$ $\displaystyle 2^{d},$ (110)
$\displaystyle\textrm{Tr}\left(h_{d,\mu}h_{d,\mu}h_{d,\nu}h_{d,\nu}\right)=$
$\displaystyle 2^{d}.$ (111)

We now consider $\textrm{Tr}\left(h_{d,\mu}h_{d,\nu}h_{d,\mu}h_{d,\nu}\right)$
with $\mu>\nu$. Let us first work out $h_{d,\mu}h_{d,\nu}$, it is given by

$\begin{split}&\sum_{\\{x\\},\\{y\\}}\overbrace{\sigma_{0}\otimes\cdots\sigma_{0}}^{d-\mu}\otimes\sigma_{\mu,\vec{x}}\otimes\overbrace{\rho_{x_{1}}\otimes\cdots\otimes\rho_{x_{\mu-\nu-1}}}^{\mu-\nu-1}\otimes\rho_{x_{\mu-\nu}}\sigma_{\nu,\vec{y}}\otimes\overbrace{\rho_{x_{\mu-\nu+1}}\rho_{y_{1}}\otimes\cdots\otimes\rho_{x_{\mu-1}}\rho_{y_{\nu-1}}}^{\nu-1}\\\
=&\sum_{\\{x\\}}\overbrace{\sigma_{0}\otimes\cdots\sigma_{0}}^{d-\mu}\otimes\sigma_{\mu,\vec{x}}\otimes\overbrace{\rho_{x_{1}}\otimes\cdots\otimes\rho_{x_{\mu-\nu-1}}}^{\mu-\nu-1}\otimes\rho_{x_{\mu-\nu}}\sigma_{\nu,\vec{y}(\vec{x})}\otimes\overbrace{\rho_{x_{\mu-\nu+1}}\otimes\cdots\otimes\rho_{x_{\mu-1}}}^{\nu-1},\end{split}$
(112)

where

$\sum_{\\{x\\}}:=\sum_{x_{1},\ldots,x_{\mu-1}},\quad\sum_{\\{y\\}}:=\sum_{y_{1},\ldots,y_{\nu-1}},$
(113)

and in the second line we have used that $y_{k}=x-\mu+k$ for
$k=1,\cdots,\nu-2$ so that

$\vec{y}(\vec{x}):=(x_{\mu-\nu+1},x_{\mu-\nu+2},\ldots,x_{\mu-1},y_{\nu}).$
(114)

Hence we can write $h_{d,\mu}h_{d,\nu}h_{d,\mu}h_{d,\nu}$ as

$\begin{split}&\sum_{\\{x\\},\\{x^{\prime}\\}}\overbrace{\sigma_{0}\otimes\cdots\sigma_{0}}^{d-\mu}\otimes\sigma_{\mu,\vec{x}}\sigma_{\mu,\vec{x}^{\prime}}\otimes\overbrace{\rho_{x_{1}}\rho_{x^{\prime}_{1}}\otimes\cdots\otimes\rho_{x_{\mu-\nu-1}}\rho_{x^{\prime}_{\mu-\nu-1}}}^{\mu-\nu-1}\\\
&\qquad\quad\otimes\rho_{x_{\mu-\nu}}\rho_{x^{\prime
c}_{\mu-\nu}}\sigma_{\nu,\vec{y}(\vec{x})}\sigma_{\nu,\vec{y}(\vec{x}^{\prime})}\otimes\overbrace{\rho_{x_{\mu-\nu+1}}\rho_{x^{\prime}_{\mu-\nu+1}}\otimes\cdots\otimes\rho_{x_{\mu-1}}\rho_{x^{\prime}_{\mu-1}}}^{\nu-1},\end{split}$
(115)

where we have used
$\sigma_{\mu,\vec{y}(\vec{x})}\rho_{x^{\prime}_{\mu-\nu}}=\rho_{x^{\prime
c}_{\mu-\nu}}\sigma_{\mu,\vec{y}(\vec{x})}$. It is clear the only nonzero
terms are those with

$x^{\prime}_{\kappa}=x_{\kappa},\quad\text{if $\kappa\neq\mu-\nu$};\quad
x^{\prime}_{\mu-\nu}=x^{c}_{\mu-\nu}=1-x_{\mu-\nu}.$ (116)

Under this condition we see (note that $\sigma_{\nu,\vec{y}(\vec{x})}$ does
not depend on the last index of $\vec{y}(\vec{x})$)

$\sigma_{\nu,\vec{y}(\vec{x})}\sigma_{\nu,\vec{y}(\vec{x}^{\prime})}=\sigma_{0},$
(117)
$\sigma_{\mu,\vec{x}}\sigma_{\mu,\vec{x}^{\prime}}=\begin{pmatrix}e^{i\phi
S_{\mu,\mu-\nu}(x_{\mu-\nu}-x_{\mu-\nu}^{c})}&0\\\ 0&e^{-i\phi
S_{\mu,\mu-\nu}(x_{\mu-\nu}-x_{\mu-\nu}^{c})}\end{pmatrix},$ (118)

whose trace is

$\textrm{Tr}\sigma_{\mu,\vec{x}}\sigma_{\mu,\vec{x}^{\prime}}=2\cos(\phi
S_{\mu,\mu-\nu})=2\cos\phi.$ (119)

Now we can perform the sum over $x$ explicitly and obtain

$\textrm{Tr}h_{d,\mu}h_{d,\nu}h_{d,\mu}h_{d,\nu}=2^{d}\cos\left(\phi
S_{\mu\nu}\right)=2^{d}\cos\phi,$ (120)

which is independent of $S_{\mu\nu}$. Combining equations (110), (111) and
(120), we obtain the total fourth moment

$\displaystyle\textrm{Tr}H^{4}=2^{d}\left[d+4\binom{d}{2}+2\binom{d}{2}\cos\phi\right]=2^{d}\left[d(d-1)(2+\cos\phi)+d\right],$
(121)

and the normalized fourth moment is equal to

$\displaystyle\frac{2^{-d}\textrm{Tr}H^{4}}{[2^{-d}\textrm{Tr}H^{2}]^{2}}=\frac{d-1}{d}(2+\cos\phi)+\frac{1}{d},$
(122)

which is in agreement with the averaged fourth moment
$\langle\textrm{Tr}H^{4}\rangle$ first obtained in Marinari:1995jwr .

###  A.2 $\textrm{Tr}H^{6}$

In this section we show that $\textrm{Tr}H^{6}$ does not depend on the
disorder realization.

Most contribution to $\textrm{Tr}H^{6}$ can be reduced to combinations
occurring in $\textrm{Tr}H^{4}$. We have two new combinations:
$\textrm{Tr}h_{\mu}h_{\nu}h_{\mu}h_{\omega}h_{\nu}h_{\omega}$ and
$\textrm{Tr}h_{\mu}h_{\nu}h_{\omega}h_{\mu}h_{\nu}h_{\omega}$ with
$\mu>\nu>\omega$. For notational clarity we focus on the cases of
$\textrm{Tr}h_{d}h_{d-1}h_{d}h_{d-2}h_{d-1}h_{d-2}$ and
$\textrm{Tr}h_{d}h_{d-1}h_{d-2}h_{d}h_{d-1}h_{d-2}$, so that we only need to
use

$\displaystyle h_{d}$ $\displaystyle=$
$\displaystyle\sum_{x_{1},\ldots,x_{d-1}}\sigma_{d,\vec{x}}\otimes\rho_{x_{1}}\otimes\rho_{x_{2}}\otimes\cdots,$
(123) $\displaystyle h_{d-1}$ $\displaystyle=$
$\displaystyle\sigma_{0}\otimes\sum_{x^{\prime}_{1},\ldots,x^{\prime}_{d-2}}\sigma_{d-1,\vec{x}^{\prime}}\otimes\rho_{x_{1}^{\prime}}\otimes\cdots,$
(124) $\displaystyle h_{d-2}$ $\displaystyle=$
$\displaystyle\sigma_{0}\otimes\sigma_{0}\otimes\sum_{x^{\prime\prime}_{1},\ldots,x^{\prime\prime}_{d-3}}\sigma_{d-2,\vec{x}^{\prime\prime}}\otimes\cdots$
(125)

with the $\cdots$ representing $d-3$ additional factors $\rho_{x}$. Each $h$
appears two times in the traces, and we use another set of dummy indices $y$,
$y^{\prime}$ and $y^{\prime\prime}$ to be the summation indices for their
second appearances. From the multiplication of the last $d-3$ factors, we know
the sum only receive contributions from

$\begin{split}x^{\prime\prime}_{1}=y^{\prime\prime}_{1}=x^{\prime}_{2}&=y^{\prime}_{2}=x_{3}=y_{3},\\\
x^{\prime\prime}_{2}=y^{\prime\prime}_{2}=x^{\prime}_{3}&=y^{\prime}_{3}=x_{4}=y_{4},\\\
&\ \vdots\\\
x^{\prime\prime}_{d-3}=y^{\prime\prime}_{d-3}=x^{\prime}_{d-2}&=y^{\prime}_{d-2}=x_{d-1}=y_{d-1},\end{split}$
(126)

and the summation symbol simplifies accordingly:

$\sum_{x_{1},\ldots,x_{d-1}}\sum_{y_{1},\ldots,y_{d-1}}\sum_{x^{\prime}_{1},\ldots,x^{\prime}_{d-2}}\sum_{y^{\prime}_{1},\ldots,y^{\prime}_{d-2}}\sum_{x^{\prime\prime}_{1},\ldots,x^{\prime\prime}_{d-3}}\sum_{y^{\prime\prime}_{1},\ldots,y^{\prime\prime}_{d-3}}\rightarrow\sum_{x_{3},\ldots,x_{d-1}}\sum_{x_{1},x_{2},x_{1}^{\prime}}\sum_{y_{1},y_{2},y_{1}^{\prime}}.$
(127)

We now work out $\textrm{Tr}(h_{d}h_{d-1}h_{d}h_{d-2}h_{d-1}h_{d-2})$. The
nontrivial part of the trace is

$\displaystyle\sum_{x_{1},x_{2},x_{1}^{\prime}}\sum_{y_{1},y_{2},y_{1}^{\prime}}\textrm{Tr}\left(\sigma_{d,\vec{x}}\sigma_{d,\vec{y}}\right)\textrm{Tr}\left(\rho_{x_{1}}\sigma_{d-1,\vec{x}^{\prime}}\rho_{y_{1}}\sigma_{d-1,\vec{y}^{\prime}}\right)\textrm{Tr}\left(\rho_{x_{2}}\rho_{x_{1}^{\prime}}\rho_{y_{2}}\sigma_{d-2,\vec{x}^{\prime\prime}}\rho_{y^{\prime}_{1}}\sigma_{d-2,\vec{y}^{\prime\prime}}\right)$
(128) $\displaystyle=$
$\displaystyle\sum_{x_{1},x_{2},x_{1}^{\prime}}\sum_{y_{1},y_{2},y_{1}^{\prime}}\textrm{Tr}\left(\sigma_{d,\vec{x}}\sigma_{d,\vec{y}}\right)\textrm{Tr}\left(\rho_{x_{1}}\rho_{y^{c}_{1}}\sigma_{d-1,\vec{x}^{\prime}}\sigma_{d-1,\vec{y}^{\prime}}\right)\textrm{Tr}\left(\rho_{x_{2}}\rho_{x_{1}^{\prime}}\rho_{y_{2}}\rho_{{y^{\prime}_{1}}^{c}}\sigma_{d-2,\vec{x}^{\prime\prime}}\sigma_{d-2,\vec{y}^{\prime\prime}}\right)$
$\displaystyle=$
$\displaystyle\sum_{x_{1},x_{2},x_{1}^{\prime}}\sum_{y_{1},y_{2},y_{1}^{\prime}}\textrm{Tr}\left(\sigma_{d,\vec{x}}\sigma_{d,\vec{y}}\right)\textrm{Tr}\left(\rho_{x_{1}}\rho_{y^{c}_{1}}\sigma_{d-1,\vec{x}^{\prime}}\sigma_{d-1,\vec{y}^{\prime}}\right)\textrm{Tr}\left(\rho_{x_{2}}\rho_{x_{1}^{\prime}}\rho_{y_{2}}\rho_{{y^{\prime}_{1}}^{c}}\right),$

where in the last equality we used
$\sigma_{d-2,\vec{x}^{\prime\prime}}\sigma_{d-2,\vec{y}^{\prime\prime}}=\sigma_{d-2,\vec{x}^{\prime\prime}}^{2}=\sigma_{0}$.
The nonzero traces are those with the extra conditions

$x_{1}=y_{1}^{c},\quad x_{2}=x_{1}^{\prime}=y_{2}={y^{\prime}_{1}}^{c}$ (129)

on top of the conditions (126). With these conditions we get

$\begin{split}\textrm{Tr}\left(\sigma_{d,\vec{x}}\sigma_{d,\vec{y}}\right)=&2\cos(\phi
S_{d,1})\\\
\sum_{x_{1}}\textrm{Tr}\left(\rho_{x_{1}}\rho_{y_{1}^{c}}\sigma_{d-1,\vec{x}^{\prime}}\sigma_{d-1,\vec{y}^{\prime}}\right)=&\textrm{Tr}\left(\sigma_{d-1,\vec{x}^{\prime}}\sigma_{d-1,\vec{y}^{\prime}}\right)=2\cos(\phi
S_{d-1,1})\\\
\textrm{Tr}\left(\rho_{x_{2}}\rho_{x_{1}^{\prime}}\rho_{y_{2}}\rho_{{y^{\prime}_{1}}^{c}}\right)=&\textrm{Tr}\left(\rho_{x_{2}}\right)=1.\end{split}$
(130)

Taking the trace over the remaining $\rho_{x_{3}}\otimes\cdots\rho_{x_{d-1}}$
and sum over the remaining indices $x_{2},x_{3},\ldots,x_{d-1}$, we finally
arrive at

$\displaystyle\textrm{Tr}h_{d}h_{d-1}h_{d}h_{d-2}h_{d-1}h_{d-2}=2^{d}\cos\left(\phi
S_{d,1}\right)\cos\left(\phi S_{d-1,1}\right)=2^{d}(\cos\phi)^{2},$ (131)

which is independent of disorder realizations.
Next we proceed to $\textrm{Tr}h_{d}h_{d-1}h_{d-2}h_{d}h_{d-1}h_{d-2}$. The
nontrivial part is given by

$\displaystyle\begin{split}&\sum_{x_{1},x_{2},x_{1}^{\prime}}\sum_{y_{1},y_{2},y_{1}^{\prime}}\textrm{Tr}\left(\sigma_{d,\vec{x}}\sigma_{d,\vec{y}}\right)\textrm{Tr}\left(\rho_{x_{1}}\sigma_{d-1,\vec{x}^{\prime}}\rho_{y_{1}}\sigma_{d-1,\vec{y}^{\prime}}\right)\textrm{Tr}\left(\rho_{x_{2}}\rho_{x_{1}^{\prime}}\sigma_{d-2,\vec{x}^{\prime\prime}}\rho_{y_{2}}\rho_{y^{\prime}_{1}}\sigma_{d-2,\vec{y}^{\prime\prime}}\right)\\\
=&\sum_{x_{1},x_{2},x_{1}^{\prime}}\sum_{y_{1},y_{2},y_{1}^{\prime}}\textrm{Tr}\left(\sigma_{d,\vec{x}}\sigma_{d,\vec{y}}\right)\textrm{Tr}\left(\rho_{x_{1}}\rho_{y^{c}_{1}}\sigma_{d-1,\vec{x}^{\prime}}\sigma_{d-1,\vec{y}^{\prime}}\right)\textrm{Tr}\left(\rho_{x_{2}}\rho_{x_{1}^{\prime}}\rho_{y_{2}^{c}}\rho_{{y^{\prime}_{1}}^{c}}\right),\end{split}$
(132)

Note only the last factor is different from that of equation (128), and the
extra conditions enforced this time are

$x_{1}=y_{1}^{c},\quad x_{2}=x_{1}^{\prime}=y_{2}^{c}={y^{\prime}_{1}}^{c}.$
(133)

So equation (132) reduces to

$\displaystyle\sum_{x_{1},x_{2}}\textrm{Tr}\left(\sigma_{d,\vec{x}}\sigma_{d,\vec{y}}\right)\textrm{Tr}\left(\rho_{x_{1}}\sigma_{d-1,\vec{x}^{\prime}}\sigma_{d-1,\vec{y}^{\prime}}\right)$
$\displaystyle=$
$\displaystyle\sum_{x_{1},x_{2}}2\cos\left[\phi(x_{1}-x_{1}^{c})S_{d,1}+\phi(x_{2}-x_{2}^{c})S_{d,2}\right]\textrm{Tr}\left[\rho_{x_{1}}\begin{pmatrix}e^{i\phi(x_{2}-x_{2}^{c})S_{d-1,1}}&0\\\
0&e^{-i\phi(x_{2}-x_{2}^{c})S_{d-1,1}}\end{pmatrix}\right]$ $\displaystyle=$
$\displaystyle 8\cos\left(\phi S_{d,1}\right)\cos\left(\phi
S_{d,2}\right)\cos\left(\phi S_{d-1,1}\right)+8i\sin\left(\phi
S_{d,1}\right)\sin\left(\phi S_{d,2}\right)\sin\left(\phi S_{d-1,1}\right).$
(134)

Note this does depend on disorder realizations of $S_{\mu\nu}$ due to the sine
terms. However, in the expansion of $\textrm{Tr}H^{6}$, the
$\textrm{Tr}h_{d}h_{d-1}h_{d-2}h_{d}h_{d-1}h_{d-2}$ can always be paired with
a reverse-ordered term, namely
$\textrm{Tr}h_{d-2}h_{d-1}h_{d}h_{d-2}h_{d-1}h_{d}$. The previous calculations
can be repeated easily, with the only change being a reverse ordering of
matrices, and instead of equation (A.2) we now have

$\displaystyle\sum_{x_{1},x_{2}}\textrm{Tr}\left(\sigma_{d,\vec{x}}\sigma_{d,\vec{y}}\right)\textrm{Tr}\left(\rho_{x_{1}}\sigma_{d-1,\vec{y}^{\prime}}\sigma_{d-1,\vec{x}^{\prime}}\right)$
$\displaystyle=$
$\displaystyle\sum_{x_{1},x_{2}}2\cos\left[\phi(x_{1}-x_{1}^{c})S_{d,1}+\phi(x_{2}-x_{2}^{c})S_{d,2}\right]\textrm{Tr}\left[\rho_{x_{1}}\begin{pmatrix}e^{-i\phi(x_{2}-x_{2}^{c})S_{d-1,1}}&0\\\
0&e^{i\phi(x_{2}-x_{2}^{c})S_{d-1,1}}\end{pmatrix}\right]$ $\displaystyle=$
$\displaystyle 8\cos\left(\phi S_{d,1}\right)\cos\left(\phi
S_{d,2}\right)\cos\left(\phi S_{d-1,1}\right)-8i\sin\left(\phi
S_{d,1}\right)\sin\left(\phi S_{d,2}\right)\sin\left(\phi S_{d-1,1}\right).$
(135)

Taking the sum of (A.2) and (A.2) we notice the sine terms cancel, thus the
result no longer depends on $S_{\mu\nu}$. After performing the sum over
$x_{3},\ldots,x_{d-1}$, we conclude

$\begin{split}&\textrm{Tr}h_{d}h_{d-1}h_{d-2}h_{d}h_{d-1}h_{d-2}+\textrm{Tr}h_{d-2}h_{d-1}h_{d}h_{d-2}h_{d-1}h_{d}\\\
=&2^{d+1}\cos\left(\phi S_{d,1}\right)\cos\left(\phi
S_{d,2}\right)\cos\left(\phi S_{d-1,1}\right)\\\
=&2^{d+1}\left(\cos\phi\right)^{3},\end{split}$ (136)

which readily generalizes to the generic cases
$\textrm{Tr}h_{\mu}h_{\nu}h_{\omega}h_{\mu}h_{\nu}h_{\omega}+\textrm{Tr}h_{\omega}h_{\nu}h_{\mu}h_{\omega}h_{\nu}h_{\mu}$
with $\mu>\nu>\omega$.

##  Appendix B Moments, words, chord diagrams and intersection graphs

In this appendix, we will discuss how the leading and subleading large $d$
contributions to moments can be obtained through chord diagrams. For the
leading contributions Parisi’s original paper Parisi:1994jg  already has a
comprehensive discussion, so we will briefly rephrase his work. Parisi also
explicitly listed the subleading contributions up to the eighteenth moment,
without giving a chord diagram interpretation of the results. We find in fact
there is a nice correspondence between subleading contributions and the
leading-contribution chord diagrams through a deletion procedure, and we will
discuss it at some length.

###  B.1 Leading contributions

The $2p$-th moment $\langle\textrm{Tr}H^{2p}\rangle$ is given by the sum of
all $2p$-step Wilson loops:

$\langle\textrm{Tr}H^{2p}\rangle=\sum_{\mathcal{C},|\mathcal{C}|=2p}\langle
W(\mathcal{C})\rangle.$ (137)

We will classify all the $2p$-step Wilson loops into groups by the total
number of Euclidean dimensions they traverse. Since the $2p$ steps need to
form a loop, at most they can traverse $p$ different dimensions. If we follow
the path of a $2p$-step loop, each time a new step is taken along a dimension
that has not been traversed, we pick up a multiplicity factor counting the
remaining dimensions. For example, the first step of any loop can freely
choose any of the $d$ dimensions; the nearest next step that takes a different
dimension has $d-1$ dimensions to choose from, and so on. By this reasoning we
see if a loop traverses $k$ dimensions, the multiplicity factor from this
effect alone is $d(d-1)\cdots(d-k+1)\sim d^{k}$. Since $k\leq p$, the leading
large $d$ contributions will come from those loops that traverse $p$ different
dimensions, having a multiplicity factor of $d(d-1)\cdots(d-p+1)$, and each of
the $p$ chosen dimensions is traversed twice to form a loop. We can use an
alphabet of $p$ different letters to represent the $p$ different dimensions,
and use a $2p$-letter word with each alphabet letter appearing twice to
represent a loop: we read the $2p$ letters in the word from left to right, and
we traverse the dimension that is represented by the letter. To avoid double
counting we should demand that the first appearances the letters in a word
must be ordered as they are in the alphabet. As a few examples, $aabb$ is a
permissible word but $bbaa$ is not; $abcacb$ is permitted but $acbabc$ is not;
$abbacc$ is permitted but $caacbb$ is not. It is easy to see there are
$(2p-1)!!$ different words we can form by having $p$ letters each appearing
twice. If we connect the same letters in a word with lines in the upper half
plane, we form what is called a chord diagram, and the lines are called the
chords. See figure 9 for a few examples.

WordLattice pathChord diagramIntersection
graph$abab$$a$$b$$a$$b$$a$$b$$a$$b$$a$$b$$abcacb$$a$$b$$c$$a$$c$$b$$a$$b$$c$$a$$c$$b$$a$$b$$c$$abcabc$$a$$b$$c$$a$$b$$c$$a$$b$$c$$a$$b$$c$$a$$b$$c$
Figure 9: Three examples of loop paths on hypercube and their representations
in terms of words, chord diagrams and intersection graphs.

The Wilson loop value can be calculated by decomposing a loop into elementary
plaquettes, namely the faces of our hypercube. We can project the path into
all the $\binom{d}{2}$ coordinate axis planes, and if the projection into the
$\mu\nu$ plane has a plaquette shape, we pick up phase of $e^{i\phi
S_{\mu\nu}}$; if the projection is a backtracking path which has zero area,
then the contribution is just $1$. Multiplying the contributions from all
projections gives the value of the Wilson loop. It is important to note for it
leading contributions the projections do not loop around the same plaquette
twice, because there are only two steps along each dimension: a step forward
and a step back. This means we cannot pick up phases like $e^{2i\phi
S_{\mu\nu}}$ from the $\mu\nu$ plane, and the disorder average over
$S_{\mu\nu}$ on each face results in a $\cos\phi$ for each projection that
loops around a face, and contributions from different plaquettes multiply.
Hence we have the following formula:

$\langle W(\mathcal{C})\rangle=q^{A(\mathcal{C})},$ (138)

with

$\begin{split}q:=&\cos\phi,\\\
A(\mathcal{C}):=&\sum_{\mu<\nu}A_{\mu\nu}(\mathcal{C}),\end{split}$ (139)

where $A_{\mu\nu}(\mathcal{C})$ is the area of the loop’s projection into
$\mu\nu$ plane which takes value of either $0$ or $1$.

In the word representation of lattice paths, if we want to study the loop
projection into a particular plane, we only need to focus on the two alphabet
letters that represent the plane. Suppose the $\mu$ and $\nu$ dimensions are
represented by the letters $a$ and $b$, respectively. To study the projection
into $\mu\nu$ plane, we can temporarily forget letters other than $a$ and $b$.
With regard to $a$ and $b$, there are only three scenarios:

$\displaystyle\ldots a\ldots a\ldots b\ldots b\ldots,$ (140)
$\displaystyle\ldots a\ldots b\ldots b\ldots a\ldots,$ (141)
$\displaystyle\ldots a\ldots b\ldots a\ldots b\ldots.$ (142)

It is clear that the first two cases have zero-area projections in the
$\mu\nu$ plane and the third case has an area-one projection. In terms of
chord diagrams, the first two have zero intersections between chords $a$ and
$b$, whereas the third has one intersection. Now we can synthesize equations
(137) and (138) for the leading contributions as

$\langle\textrm{Tr}H^{2p}\rangle_{\text{leading}}=d(d-1)\cdots(d-p+1)\sum_{C(p\text{-chord})}q^{\\#\text{
of intersections in }C},$ (143)

where $q=\cos\phi$ and $C$ denotes chord diagrams with $p$ chords. In other
words the leading moments are the generating functions of chord intersections.
The moments calculated by (143) also appear in the Sachdev-Ye-Kitaev model
Erdos:2014a ; Garcia-Garcia:2017pzl ; Cotler2016 ; Garcia-Garcia:2018kvh . The
sum on the right-hand side of (143) has an interesting solution: in his
original paper Parisi Parisi:1994jg  already suggested mapping the sum to the
vacuum expectation values of some observables in the $q$-deformed harmonic
oscillator system. This approach was further elaborated in Marinari:1995jwr .
In fact, much earlier this chord diagram sum was studied by Touchard
Touchard:1952a  and Riordan Riordan:1975a  in a more combinatorial vein, which
led to the Riordan-Touchard formula:

$\sum_{C(p\text{-chord})}q^{\\#\text{ of intersections in
}C}=\frac{1}{(1-q)^{p}}\sum\limits_{k=-p}^{p}(-1)^{k}q^{k(k-1)/2}\binom{2p}{p+k}.$
(144)

Using this formula we can quite effortlessly generate

$\displaystyle\langle\textrm{Tr}H^{2}\rangle_{\text{leading}}=$ $\displaystyle
d,$ (145) $\displaystyle\langle\textrm{Tr}H^{4}\rangle_{\text{leading}}=$
$\displaystyle d(d-1)\left[2+q\right],$ (146)
$\displaystyle\langle\textrm{Tr}H^{6}\rangle_{\text{leading}}=$ $\displaystyle
d(d-1)(d-2)\left[5+6q+3q^{2}+q^{3}\right],$ (147)
$\displaystyle\langle\textrm{Tr}H^{8}\rangle_{\text{leading}}=$ $\displaystyle
d(d-1)(d-2)(d-3)\left[14+28q+28q^{2}+20q^{3}+10q^{4}+4q^{5}+q^{6}\right],$
(148)

and so on.

###  B.2 Subleading contributions

The goal of this section is to give a chord diagram interpretation of the
subleading contributions to moments. To be clear, there are already
contributions subleading in $d$ included in equation (143) due to the
multiplicity factor $d(d-1)\cdots(d-p+1)$. However, there are still subleading
contributions from $2p$-step loops that traverse only $p-1$ dimensions which
gives a multiplicity factor $d(d-1)\cdots(d-p+2)$. This section will be about
such loops. We first demonstrate that there is a bijection between subleading
words and certain structures of the leading words. The choice for this
bijection is not unique, different choices lead to different schemes of
calculating the subleading contributions, and unsurprisingly all schemes give
the same result.

####  B.2.1 The interlace scheme

As already discussed, the leading words with $2p$ letters are words with $p$
different pairs of alphabet letters. By the previous discussion it is clear
the subleading words with $2p$ letters have one alphabet letter appearing four
times, and $p-2$ other alphabet letters each appearing twice. This reflects
the fact that subleading Wilson loops discussed at the beginning of last
section must traverse a dimension four times and other remaining dimensions
two times each. For examples $aaaabb$, $abbbab$ and $abbabb$ are some
subleading words for $p=3$. From the general formula777We can write down the
general formula after some thought. The total number of $2p$-letter words
formed by $k$ alphabet letters (each can appear even number of times) is
$\sum_{[m_{1},m_{2},\ldots,m_{k}]_{p}}\binom{2p-1}{2m_{1}-1}\binom{2p-2m_{1}-1}{2m_{2}-1}\cdots\binom{2p-2\sum_{l=1}^{k-1}m_{l}-1}{2m_{k}-1},$
(149) where $[m_{1},m_{2},\ldots,m_{k}]_{p}$ denotes a composition of $p$,
that is, an ordered $k$-tuple $(m_{1},m_{2},\ldots,m_{k})$ such that
$\sum_{l=1}^{k}m_{l}=p$. it is clear that we can form
$\frac{1}{3}\binom{p}{2}(2p-1)!!$ subleading words of length $2p$.

The following map is a bijection between subleading words and interlacing
structures of leading words:

$\ldots a\ldots a\ldots a\ldots
a\ldots\mapsto\ldots\underline{a}\ldots\underline{b}\ldots\underline{a}\ldots\underline{b}\ldots,$
(150)

where the $\cdots$ part remains unchanged after the mapping. We added
underlines on the right to emphasize the map is toward an interlacing
structure, instead of the leading word that contains this interlacing
structure. In the context of this mapping, it is convenient for us to adopt a
“jump an alphabet letter” convention for subleading words: we jump over the
alphabet letter that immediately follows (in the alphabet) the letter that
appears four times in the word. For example, $aaaabb$ and $aaaacc$ are
equivalent words, but we prefer the second representation because it is mapped
to $ababcc$ without changing the letter $c$. Let us also see an example of the
inverse mapping. The leading word $abcacb$ has two interlacing structures,
each will be mapped to a subleading word:

$\displaystyle\underline{a}\underline{b}c\underline{a}c\underline{b}$
$\displaystyle\mapsto aacaca,$ (151)
$\displaystyle\underline{a}b\underline{cac}b$ $\displaystyle\mapsto abaaab.$
(152)

It is clear the mappings to both directions are injective and hence bijective,
because two different objects before mapping can only be different by virtue
of the $\cdots$ parts in Eq. 150, the difference is inherited after the
mapping. Note each interlacing structure in a leading word corresponds to an
intersection in the corresponding chord diagram, so we may also say there is a
bijection between subleading words and the intersections of the leading chord
diagrams.888A byproduct of this discussion is that we just completed a
bijective proof of the following statement: the total number of intersections
among all chord diagrams with $p$ chord is $\frac{1}{3}\binom{p}{2}(2p-1)!!$.
We can easily generalize the proof to other intersection structures. Other
proofs of this statement already exist, see for example Flajolet:1997a ;
Jia:2018ccl .

We have demonstrated that the bijection (150) allows us to use the interlacing
structures in leading words to represent the subleading Wilson loops. The
remaining question is how to read off the values of the Wilson loops from the
leading word interlacing structures. Let us recall that for a leading Wilson
loop, each interlacing structure in its word representation represents a
projection of the path that loops around a plaquette. Obviously, after the
mapping (150), this particular interlacing structure is removed, and the
leading Wilson loop becomes a subleading Wilson loop in which this plaquette
projection gets squashed to a zero-area projection. However, this is not the
end of the story: it is conceivable that the removal of one interlacing
structure interferes with other interlacing structures in the same word, so
that more plaquette-shaped projections get squashed as a result. We are faced
with three possibilities:

  1. 1.

The other interlacing structure that might be interfered with by the removal
of $\underline{a}\underline{b}\underline{a}\underline{b}$ is formed by two
other alphabet letters, so we have

$\underline{a}\ldots\underline{b}\ldots\underline{a}\ldots\underline{b}\ldots
c\ldots d\ldots c\ldots d\ldots\mapsto a\ldots a\ldots a\ldots a\ldots c\ldots
d\ldots c\ldots d\ldots.$ (153)

In this case the plaquette projection represented by $cdcd$ cannot be affected
because the hypercube dimensions represented by $c$ and $d$ are in the
orthogonal complement of $a$ and $b$.

  2. 2.

The other interlacing structure that might be interfered with by the removal
of $\underline{a}\underline{b}\underline{a}\underline{b}$ is formed by one
other alphabet letter interlacing with one of $a$ or $b$ but not both, so we
have

$\ldots\underline{a}\ldots\underline{b}\ldots c\ldots\underline{a}\ldots
c\ldots\underline{b}\ldots\mapsto\ldots a\ldots a\ldots c\ldots a\ldots
c\ldots a\ldots.$ (154)

The effect on the corresponding Wilson loops can be read off visually:

$a$$b$$c$$a$$c$$b$$\mapsto$$a$$a$$a$$c$$a$$c$ (155)

and we see the projected area in the $ac$ plane remains one.

  3. 3.

Another alphabet letter $c$ interlaces with both $a$ and $b$, and we want to
investigate what happens to the projected areas represented by the two new
interlacing structures $acac$ and $bcbc$. The word map is

$\ldots\underline{a}\ldots\underline{b}\ldots
c\ldots\underline{a}\ldots\underline{b}\ldots c\ldots\mapsto\ldots a\ldots
a\ldots c\ldots a\ldots a\ldots c\ldots.$ (156)

The corresponding Wilson loop transforms as

$a$$b$$c$$a$$b$$c$$\mapsto$$a$$a$$a$$a$$c$$c$ (157)

and we see the all three plaquettes in the Wilson loop before the mapping
collapse to zero area after the mapping.

We can summarize the above three cases as the following: for any subleading
Wilson loop represented by the leading word interlacing structure
$\ldots\underline{a}\ldots\underline{b}\ldots\underline{a}\ldots\underline{b}\ldots$,
this subleading Wilson loop has the value

$q^{\text{\\# of interlacing structures in this word}-2(\text{\\# of
triangular structures containing $\underline{abab}$)}-1},$ (158)

where “triangular structures containing $\underline{abab}$” are structures
like $\ldots\underline{a}\ldots\underline{b}\ldots
c\ldots\underline{a}\ldots\underline{b}\ldots c\ldots$. For example, if the
leading word is $abcdabcd$ and the interlacing structure we are interested in
$\underline{ab}cd\underline{ab}cd$, then there are two triangular structures
containing $\underline{abab}$, namely
$\underline{ab}c\cdot\underline{ab}c\cdot$ and $\underline{ab}\cdot
d\underline{ab}\cdot d$. There are six interlacing structures in the word
$abcdabcd$, so the subleading Wilson loop has the value $q^{6-(2\times
2+1)}=q$.

We can obtain a rather compact and visual representation of rule (158) if we
introduce the notion of intersection graphs. To obtain the intersection graph
of a leading word, we first draw its chord diagram. The intersection graph is
then obtained by the following two steps:

  1. 1.

represent every chord by a vertex,

  2. 2.

connect two vertices if and only if the chords they represent intersect each
other.

Intersection graph |  |  |  |
---|---|---|---|---
Leading value | 1 | $q$ | $q^{2}$ | $q^{3}$
Subleading value | 0 | $1$ | $2q$ | $3$
Multiplicity | 5 | 6 | 3 | 1
Table 2: All the intersection graphs for the sixth moment. Intersection graph |  |  |  |  |  |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---|---|---
Leading value | 1 | $q$ | $q^{2}$ | $q^{2}$ | $q^{3}$ | $q^{3}$ | $q^{3}$ | $q^{4}$ | $q^{4}$ | $q^{5}$ | $q^{6}$
Subleading value | 0 | $1$ | $2q$ | $2q$ | $3q^{2}$ | $3q^{2}$ | $3$ | $q^{3}+3q$ | $4q^{3}$ | $4q^{2}+1$ | $6q$
Multiplicity | 14 | 28 | 4 | 24 | 4 | 8 | 8 | 8 | 2 | 4 | 1
Table 3: All the intersection graphs for the eighth moment.

We refer readers to figure 9 for a few examples. In the intersection graph
language, the leading moment (143) can be written as

$\langle\textrm{Tr}H^{2p}\rangle_{\text{leading}}=d(d-1)\cdots(d-p+1)\sum_{G(p\text{-vertex})}q^{E_{G}},$
(159)

where the sum is over all the $(2p-1)!!$ intersection graphs $G$, and $E_{G}$
denotes the total number of edges in $G$. And from formula (158), the
subleading moment can be written as

$\langle\textrm{Tr}H^{2p}\rangle_{\text{subleading}}=d(d-1)\cdots(d-p+2)\sum_{G(p\text{-vertex})}\sum_{e\in
G}q^{E_{G}-(2T_{e}+1)},$ (160)

where $e$ denotes edges in $G$ and $T_{e}$ is the number of triangles that has
$e$ as one of its sides. Notice in intersection graphs, the triangular
structures in words literally become triangles. So equation (160) is telling
us to go through all the edges of the intersection graphs one by one, delete
the edge we are looking at and all the triangles that has it as a side, then
count the number of edges of the remaining graph, and that is the power we
raise $q$ to. Let us work out how equation (160) for low-order moments: in
tables 2 and 3, all the intersection graphs contributing to the sixth and the
eighth moment are respectively listed. The leading-contribution values they
represent are just $q$ raised to the powers being the numbers of edges of
those graphs. The subleading-contribution values are obtained by the edge and
triangle deletion procedure just described. After summing over all graphs we
can check the total leading contributions are just those given by equations
(147) and (148); the subleading contributions are

$\displaystyle\langle\textrm{Tr}H^{6}\rangle_{\text{subleading}}=$
$\displaystyle d(d-1)\left[9+6q\right],$ (161)
$\displaystyle\langle\textrm{Tr}H^{8}\rangle_{\text{subleading}}=$
$\displaystyle d(d-1)(d-2)\left[56+86q+52q^{2}+16q^{3}\right],$ (162)

which are consistent with the results of Marinari, Parisi and
RitortMarinari:1995jwr . For subleading contributions of higher moments, we
refer readers to the same reference.

It would be very useful to develop a Riordan-Touchard-like formula for
subleading moments (160), but we have not found one yet.

####  B.2.2 The nest scheme and the alignment scheme

The readers may have noticed that we can easily form two other bijections
similar to equation (150), namely:

$\ldots a\ldots a\ldots a\ldots
a\ldots\mapsto\ldots\underline{a}\ldots\underline{b}\ldots\underline{b}\ldots\underline{a}\ldots$
(163)

or

$\ldots a\ldots a\ldots a\ldots
a\ldots\mapsto\ldots\underline{a}\ldots\underline{a}\ldots\underline{b}\ldots\underline{b}\ldots.$
(164)

In some literature kim:2015a  the $abba$ structure is called a nest and the
$aabb$ structure is called an alignment. Hence we will call the calculations
based on the former the nest scheme and the latter the alignment scheme. By
the same reasoning in the interlace scheme section, we know there are exactly
the same total number of interlaces, nests and alignments when all the chord
diagrams with $p$ chords are counted, which is
$\frac{1}{3}\binom{p}{2}(2p-1)!!$. We can do a hypercube Wilson loop analysis
similar to that of the interlace scheme, namely the analysis wrapping around
equations (153)-(156) and see what the reduction to subleading words does to
the power of $q$. The end result is the following: for both the nest scheme
and the alignment scheme when there is a third chord (in terms of chord
diagrams) intersecting both chords represented by
$\underline{a}\ldots\underline{b}\ldots\underline{b}\ldots\underline{a}$ (nest
scheme) or
$\underline{a}\ldots\underline{a}\ldots\underline{b}\ldots\underline{b}$
(alignment scheme), the power of $q$ reduces by two. In all other scenarios
the power of $q$ remains the same. That is,

  1. 1.

The nest scheme:

$\begin{split}\ldots\underline{a}\ldots\underline{b}\ldots
c\ldots\underline{b}\ldots\underline{a}\ldots
c&\mapsto\ldots\underline{a}\ldots\underline{a}\ldots
c\ldots\underline{a}\ldots\underline{a}\ldots c\ldots\\\ \implies
q^{\\#}&\mapsto\ q^{\\#-2},\end{split}$ (165)

otherwise $q^{\\#}\mapsto\ q^{\\#}.$

  2. 2.

The alignment scheme999Here we temporarily use $\underline{adad}$ instead of
$\underline{abab}$, so that upon the insertion of $c$, we still comply with
the convention that what comes earlier in the alphabet comes earlier in the
word.

$\begin{split}\ldots\underline{a}\ldots
c\ldots\underline{a}\ldots\underline{d}\ldots
c\ldots\underline{d}\ldots&\mapsto\ldots\underline{a}\ldots
c\ldots\underline{a}\ldots\underline{a}\ldots c\ldots\underline{a}\ldots\\\
\implies q^{\\#}&\mapsto\ q^{\\#-2},\end{split}$ (166)

otherwise $q^{\\#}\mapsto\ q^{\\#}.$

In terms intersection graphs, we cannot distinguish a nest from an alignment
because both are represented by a pair of vertices not connected by any edge.
Hence in terms of intersection graphs we can at best give a prescription in
terms of the sum of the nest scheme and the alignment scheme, which gives two
times the subleading contribution. A chord that intersects both chords of a
nest or an alignment translates to a “wedge” structure in intersections
graphs, see figure 10. Therefore, the prescription for the sum of the nest
scheme and alignment scheme is this: for all the pairs of the vertices that
are not connected by any edge in an intersection graph, delete all the
“wedges” that connect the two vertices. The number of edges in the resulting
graph is the power on $q$. The sum of all such resulting graphs from all
leading intersection graphs gives two times the subleading coefficients of
moments.

$\underline{a}$$\underline{b}$$c$$\underline{b}$$\underline{a}$$c$nest

$\underline{a}$$c$$\underline{a}$$\underline{d}$$c$$\underline{d}$alignment

$c$$\underline{a}$$\underline{b}$ (or $\underline{d})$ wedge

Figure 10: The two left figures show a chord $c$ intersecting both chords in a
nest structure and and in an alignment structure. Note only the underlined
letters represent the nest or alignment structures. The right figure shows the
intersection graph of such scenarios.

####  B.2.3 The averaged scheme

The interlace scheme picks all the edges in the leading intersection graphs,
whereas the nest and the alignment schemes pick all the pairs of vertices not
connected by any edge. All three schemes give the same contribution, so we can
average over all three schemes and get a prescription that picks all pairs of
vertices in intersection graphs, regardless of whether the pairs are connected
by any edge or not. It is clear we can combine all the scheme prescriptions
into the following one: for every pair of vertices in a leading intersection
graph, delete the edge that connects the two vertices if there is one, and
delete all the wedges that has the two vertices as the two ends.101010For
example $\underline{a}$ and $\underline{b}$ are two ends of the wedge in the
third figure in figure 10, whereas $c$ is not an end of the wedge. Raise $q$
to power of the number of edges of the resulting graph and sum them over all
such graphs. The subleading coefficient is one third of this sum.

There is a more graphical way to describe the edge and wedge deletion
prescription. That is, we take a pair of vertices and merge them into one
vertex, and all the edges before merging are inherited. However, loops (edge
that connects a vertex to itself) and double edges (two edges connecting the
same two vertices) may appear after merging, and we delete all the loops and
double edges to form a subleading intersection graph. Note that deleting a
loop is equivalent to deleting the edge connecting a chosen pair in the
language of the last paragraph, and deleting a double edge is equivalent to
deleting the wedge. Figure 11 demonstrates a few examples of such “merge and
delete” process. We can summarize the averaged scheme into one formula:

$\frac{\langle\textrm{Tr}H^{2p}\rangle_{\text{subleading}}}{d(d-1)\cdots(d-p+2)}=\frac{1}{3}\sum_{G}\sum_{\\{v_{1},v_{2}\\}\subset
v(G)}q^{E\left[G_{(v_{1},v_{2})}\right]},$ (167)

where $G$’s are all the intersection graphs formed by all the chord diagrams
with $p$ chords; $v_{1},v_{2}$ are any two vertices of $G$ and $v(G)$ denotes
the vertex set of $G$; $G_{(v_{1},v_{2})}$ is the graph formed by the “merge
and delete” procedure applied to $G$ with respect to $v_{1}$ and $v_{2}$
(namely, $v_{1}$ and $v_{2}$ are merged), and finally
$E\left[G_{(v_{1},v_{2})}\right]$ is the number of edges in
$G_{(v_{1},v_{2})}$.

$a$$b$$c$$\rightarrow$$a$$c$$\rightarrow$$a$$c$merge $ab:$

$a$$b$$c$$\rightarrow$$a$$c$$\rightarrow$$a$$c$merge $bc:$

$a$$b$$c$$\rightarrow$$a$$c$$\rightarrow$$a$$c$merge $bc:$

Figure 11: Three examples of the “merge and delete” prescription. In the first
figure, a loop is formed after merging and then deleted, in the end $q^{2}$ is
reduced to $q$; in the second figure, a double edge is formed after merging
and then deleted, and $q^{2}$ is reduced to $1$; in the third figure, a loop
and a double edge are formed and deleted, as a result $q^{3}$ is reduced to
$1$.

## References

  * [1] Subir Sachdev and Jinwu Ye.  Gapless spin-fluid ground state in a random quantum heisenberg magnet.  Phys. Rev. Lett., 70:3339–3342, May 1993.
  * [2] Alexander Kitaev.  A simple model of quantum holography.  KITP strings seminar and Entanglement 2015 program, 12 February, 7 April and 27 May 2015, http://online.kitp.ucsb.edu/online/entangled15/.
  * [3] K.K Mon and J.B French.  Statistical properties of many-particle spectra.  Annals of Physics, 95(1):90 – 111, 1975.
  * [4] T. A. Brody, J. Flores, J. B. French, P. A. Mello, A. Pandey, and S. S. M. Wong.  Random-matrix physics: spectrum and strength fluctuations.  Rev. Mod. Phys., 53:385–479, Jul 1981.
  * [5] L Benet and H A Weidenmüller.  Review of the k -body embedded ensembles of gaussian random matrices.  Journal of Physics A: Mathematical and General, 36(12):3569, 2003\.
  * [6] F. Borgonovi, F.M. Izrailev, L.F. Santos, and V.G. Zelevinsky.  Quantum chaos and thermalization in isolated systems of interacting particles.  Physics Reports, 626:1 – 58, 2016.  Quantum chaos and thermalization in isolated systems of interacting particles.
  * [7] Fausto Borgonovi, Felix M. Izrailev, and Lea F. Santos.  Timescales in the quench dynamics of many-body quantum systems: Participation ratio versus out-of-time ordered correlator.  Phys. Rev., E99(5):052143, 2019.
  * [8] A. Georges, O. Parcollet, and S. Sachdev.  Quantum fluctuations of a nearly critical heisenberg spin glass.  Physical Review B, 63(13), Mar 2001.
  * [9] Subir Sachdev.  Bekenstein-hawking entropy and strange metals.  Phys. Rev. X, 5:041025, Nov 2015.
  * [10] Juan Maldacena and Douglas Stanford.  Comments on the sachdev-ye-kitaev model.  arXiv preprint arXiv:1604.07818, 2016.
  * [11] Juan Maldacena, Stephen H. Shenker, and Douglas Stanford.  A bound on chaos.  JHEP, 08:106, 2016.
  * [12] O. Bohigas, M. J. Giannoni, and C. Schmit.  Characterization of chaotic quantum spectra and universality of level fluctuation laws.  Phys. Rev. Lett., 52:1–4, Jan 1984.
  * [13] T.H. Seligman, J.J.M. Verbaarschot, and M.R. Zirnbauer.  Quantum spectra and transition from regular to chaotic classical motion.  Physical Review Letters, 53, 1984.
  * [14] Antonio M. García-García and Jacobus J. M. Verbaarschot.  Spectral and thermodynamic properties of the Sachdev-Ye-Kitaev model.  Phys. Rev., D94(12):126010, 2016.
  * [15] Jordan S. Cotler, Guy Gur-Ari, Masanori Hanada, Joseph Polchinski, Phil Saad, Stephen H. Shenker, Douglas Stanford, Alexandre Streicher, and Masaki Tezuka.  Black holes and random matrices.  Journal of High Energy Physics, 05(5):118, 2017.
  * [16] Phil Saad, Stephen H. Shenker, and Douglas Stanford.  A semiclassical ramp in SYK and in gravity.  2018\.
  * [17] Alexander Altland and Dmitry Bagrets.  Quantum ergodicity in the SYK model.  Nucl. Phys., B930:45–68, 2018.
  * [18] Yiyang Jia and Jacobus J. M. Verbaarschot.  Spectral Fluctuations in the Sachdev-Ye-Kitaev Model.  2019\.
  * [19] J. Flores, M. Horoi, M. Müller, and T. H. Seligman.  Spectral statistics of the two-body random ensemble revisited.  Physical Review E, 63(2), Jan 2001.
  * [20] Antonio M. García-García, Yiyang Jia, and Jacobus J. M. Verbaarschot.  Universality and Thouless energy in the supersymmetric Sachdev-Ye-Kitaev Model.  Phys. Rev., D97(10):106003, 2018.
  * [21] Hrant Gharibyan, Masanori Hanada, Stephen H. Shenker, and Masaki Tezuka.  Onset of Random Matrix Behavior in Scrambling Systems.  JHEP, 07:124, 2018.  [Erratum: JHEP02,197(2019)].
  * [22] Edward Witten.  An SYK-Like Model Without Disorder.  J. Phys., A52(47):474002, 2019.
  * [23] Igor R. Klebanov and Grigory Tarnopolsky.  Uncolored random tensors, melon diagrams, and the Sachdev-Ye-Kitaev models.  Phys. Rev., D95(4):046004, 2017.
  * [24] Igor R. Klebanov, Preethi N. Pallegar, and Fedor K. Popov.  Majorana Fermion Quantum Mechanics for Higher Rank Tensors.  Phys. Rev., D100(8):086003, 2019.
  * [25] Jaewon Kim, Igor R. Klebanov, Grigory Tarnopolsky, and Wenli Zhao.  Symmetry Breaking in Coupled SYK or Tensor Models.  Phys. Rev., X9(2):021043, 2019.
  * [26] Chethan Krishnan, K. V. Pavan Kumar, and Dario Rosa.  Contrasting SYK-like Models.  JHEP, 01:064, 2018.
  * [27] Giorgio Parisi.  D-dimensional arrays of Josephson junctions, spin glasses and q deformed harmonic oscillators.  J. Phys. A: Math. Gen., 27:7555–7568, 1994.
  * [28] Enzo Marinari, Giorgio Parisi, and Felix Ritort.  Replica Theory and Large D Josephson Junction Hypercubic Models.  J. Phys. A: Math. Gen., 28:4481–4503, 1995.
  * [29] Andrea Cappelli and Filippo Colomo.  Solving the frustrated spherical model with q polynomials.  J. Phys. A, 31:3141–3151, 1998.
  * [30] Filippo Colomo.  Counting non-planar diagrams: an exact formula.  Physics Letters A, 284(1):12 – 15, 2001.
  * [31] Filippo Colomo.  Area versus length distribution for closed random walks.  J. Phys. A, 36:1539–1552, 2003.
  * [32] Antonio M. García-García, Bruno Loureiro, Aurelio Romero-Bermúdez, and Masaki Tezuka.  Chaotic-Integrable Transition in the Sachdev-Ye-Kitaev Model.  Phys. Rev. Lett., 120(24):241603, 2018.
  * [33] Tomoki Nosaka, Dario Rosa, and Junggi Yoon.  The Thouless time for mass-deformed SYK.  JHEP, 09:041, 2018.
  * [34] Tomoki Nosaka and Tokiro Numasawa.  Quantum Chaos, Thermodynamics and Black Hole Microstates in the mass deformed SYK model.  12 2019.
  * [35] Juan Maldacena and Xiao-Liang Qi.  Eternal traversable wormhole.  2018\.
  * [36] J. Zak.  Magnetic Translation Group.  Phys. Rev., 134:A1602–A16, 1964.
  * [37] R. Rammal and J. Bellissard.  An algebraic semi-classical approach to Bloch electrons in a magnetic field.  J. Phys. France, 51:1803–1830, 1990.
  * [38] P. B. Wiegmann and A. V. Zabrodin.  Quantum group and magnetic translations. Bethe Ansatz for Asbel-Hofstadter problem.  Nucl. Phys., B422:495–514, 1994.
  * [39] Dmitry Bagrets, Alexander Altland, and Alex Kamenev.  Sachdev–ye–kitaev model as liouville quantum mechanics.  Nuclear Physics B, 911:191–205, 2016.
  * [40] Antonio M. García-García and Jacobus J. M. Verbaarschot.  Analytical Spectral Density of the Sachdev-Ye-Kitaev Model at finite N.  Phys. Rev., D96(6):066012, 2017.
  * [41] Yiyang Jia and Jacobus J. M. Verbaarschot.  Large $N$ expansion of the moments and free energy of Sachdev-Ye-Kitaev model, and the enumeration of intersection graphs.  JHEP, 11:031, 2018.
  * [42] Micha Berkooz, Prithvi Narayan, and Joan Simon.  Chord diagrams, exact correlators in spin glasses and black hole bulk reconstruction.  JHEP, 08:192, 2018.
  * [43] Micha Berkooz, Mikhail Isachenkov, Vladimir Narovlansky, and Genis Torrents.  Towards a full solution of the large N double-scaled SYK model.  JHEP, 03:079, 2019.
  * [44] Wenbo Fu, Davide Gaiotto, Juan Maldacena, and Subir Sachdev.  Supersymmetric Sachdev-Ye-Kitaev models.  Phys. Rev. D, 95(2):026009, 2017.  [Addendum: Phys.Rev.D 95, 069904 (2017)].
  * [45] Takuya Kanazawa and Tilo Wettig.  Complete random matrix classification of SYK models with $\mathcal{N}=0$, $1$ and $2$ supersymmetry.  JHEP, 09:050, 2017.
  * [46] Antonio M. García-García, Tomoki Nosaka, Dario Rosa, and Jacobus J. M. Verbaarschot.  Quantum chaos transition in a two-site Sachdev-Ye-Kitaev model dual to an eternal traversable wormhole.  Phys. Rev., D100(2):026002, 2019.
  * [47] Ken-ichi Sekiguchi, Tomohiro Okamoto, and Takanori Fujiwara.  Magnetic translation symmetry on the lattice.  Prog. Theor. Phys., 122:293–307, 2009.
  * [48] Erdös, L. and Schröder, D.  Phase Transition in the Density of States of Quantum Spin Glasses.  Math. Phys. Anal. Geom., pages 441–464, 2014.
  * [49] M.E.H. Ismail, Stanton C., and Viennot G.  The Combinatorics of q-Hermite polynomials and the Askey-Wilson Integral.  Europ. J. Combinatorics, 8:379, 1987.
  * [50] A. Delon, R. Jost, and M. Lombardi.  NO2 jet cooled visible excitation spectrum - vibronic chaos induced by the x2a1-a2b2 interaction.  J. Chem. Phys., 95:5701–5718, 1991.
  * [51] Jacobus J. M. Verbaarschot.  The Spectrum of the QCD Dirac operator and chiral random matrix theory: The Threefold way.  Phys. Rev. Lett., 72:2531–2533, 1994.
  * [52] Edward V. Shuryak and J.J.M. Verbaarschot.  Random matrix theory and spectral sum rules for the Dirac operator in QCD.  Nucl. Phys. A, 560:306–320, 1993.
  * [53] J.J.M. Verbaarschot and T. Wettig.  Random matrix theory and chiral symmetry in QCD.  Ann. Rev. Nucl. Part. Sci., 50:343–410, 2000.
  * [54] J.J.M. Verbaarschot and I. Zahed.  Spectral density of the QCD Dirac operator near zero virtuality.  Phys. Rev. Lett., 70:3852–3855, 1993.
  * [55] D. Toublan and J.J.M. Verbaarschot.  Statistical properties of the spectrum of the QCD Dirac operator at low-energy.  Nucl. Phys. B, 603:343–368, 2001.
  * [56] William Cottrell, Ben Freivogel, Diego M. Hofman, and Sagar F. Lokhande.  How to Build the Thermofield Double State.  JHEP, 02:058, 2019.
  * [57] Juan Martin Maldacena.  Eternal black holes in anti-de Sitter.  JHEP, 04:021, 2003.
  * [58] A. del Campo, J. Molina-Vilaplana, L. F. Santos, and J. Sonner.  Decay of a Thermofield-Double State in Chaotic Quantum Systems.  Eur. Phys. J. ST, 227(3-4):247–258, 2018.
  * [59] Antonio M. García-García, Yiyang Jia, and Jacobus J. M. Verbaarschot.  Exact moments of the Sachdev-Ye-Kitaev model up to order $1/N^{2}$.  JHEP, 04:146, 2018.
  * [60] J. Touchard.  Sur un probleme de configurations et sur les fractions continues.  Canadian J. Math., 4:2–25, 1952.
  * [61] J. Riordan.  The distribution of crossings of chords joining pairs of $2n$ points on a circle.  Math. Comp., 29:215–222, 1975.
  * [62] Philippe Flajolet and Marc Noy.  Analytic Combinatorics of Non-crossing Configurations.  Research Report, RR-3196, INRIA, 1997.
  * [63] Jang Soo Kim and Dennis Stanton.  Bootstrapping and askey-wilson polynomials.  Journal of Mathematical Analysis and Applications, 421:501–520, 2014.