Chaos on the hypercube

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 the 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.


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 [1,2] formerly known as the two-body random ensemble [3], has greatly improved our understanding of the relationship between many-body chaos, disorder and spectral properties of the underlying Hamiltonian (see [4][5][6][7] 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 [8]. For the same reason, the SYK model can be used to address questions related to understanding micro-states and entropy of black holes [9]. There are different ways to measure the chaotic properties of the SYK model. The shorttime behavior of the Out-of-Time-Order Correlator (OTOC), which in the classical limit describes the exponential divergence of classical trajectories, was shown [10] to saturate the chaos bound [11]. 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 [12,13]. Indeed this was confirmed by numerical and analytical studies of the SYK model [14][15][16][17][18]. 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 [19] 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 [16-18, 20, 21] 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. 1 This gives an estimate of 2 N/2 /N for the spectral range of RMT fluctuations or a times scale of N 2 −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 [22] which also has been confirmed by melonic models which have similar properties in the absence of disorder [23][24][25][26]. 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 [27][28][29][30][31] 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 [32][33][34] 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 φ. At φ = 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 [35]. The spectral flow as a function of φ 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 [36][37][38], 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 φ. 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 [15,[39][40][41][42][43]. For φ > π/2, Q becomes negative and spectrum splits into two bands, which also happens for the supercharge of the supersymmetric SYK model [20,44,45]. 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 [16][17][18][19][20][21]. 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 φ < π/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 [27] 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 x with components x µ ∈ {0, 1}. The model considers a constant magnetic field such that the fluxes through all faces of the hypercube have the same magnitude φ, but with random orientations. That is, we have the field strength tensor where S µν is an antisymmetric tensor with random entries ±1 with equal probabilities. Hence we are dealing with a finite ensemble with 2 ( d 2 ) disorder realizations. We can work in the axial gauge so that the link variables are given by which is the phase we associate with the link emanating from site x along the µ-th direction. Note the sum is over all the ν's with ν < µ, and if µ = 1 we define U 1 ( 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 whereê µ is the unit basis vector in the µ-th direction. When φ = 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 where ϕ is a scalar quantum field. However, in this paper we take a first quantized view and concern ourselves with the H defined in equation (2.3), and the wave functions live in C 2 d . Let us be very explicit on how to write the Hamiltonian matrix as a two-dimensional array of numbers: since x = (x 1 , x 2 , . . . , x d ) is a string of 0 and 1's of length d, we can naturally think of 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 ∈ {1, 2, 3, . . . , 2 d } through the relation where [m − 1] 2 denotes the number m − 1 in the binary representation, and . . 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 , . . . , 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 and so on.
The Hamiltonian can be obtained recursively: where C d−1 is a diagonal unitary matrix with entries where x ν (k) is the ν-th digit of [k − 1] 2 , as defined in equation (2.5). We can verify that the following relation holds: which will be useful for section 4. For later convenience, we also introduce the notation S ρν , S 1 := 0, (2.10) so that the right-hand side of equation (2.9) is simply e iφS d .

Tensor product representation of the Hamiltonian
Since the interaction between two sites can be written in terms of the Pauli matrix σ 1 , it is not surprising that the Hamiltonian can be expressed in terms of tensor products of Pauli σ matrices. For d = 2 it is clear from equation (2.6) that Notice that the definition of σ d, x does not depend on the last component x d of x, for example we have For higher dimensions we have and in general we have . (3.6)

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 γ matrices, which then can be interpreted and the Hamiltonian for a system of 2d Majorana fermions. The simples case is φ = 0. Then the Hamiltonian is just the adjacency matrix of a hypercube graph. In the tensor product representation it is given by If we introduce the gamma matrices then the Hamiltonian can be written as This is exactly the interaction term in the Maldacena-Qi model [35]. This interaction term was shown [46] to have the spectrum Indeed this is also the well-known spectrum for the hypercube adjacency matrix. At φ = 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 γ matrices appear in the Hamiltonian, which makes the Hamiltonian look very complicated.

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 (2.3) has a sublattice symmetry so that all eigenvalues appear in pairs ±λ k . In the tensor product representation described in section 3, Γ 5 has the simple form of a tensor product of σ 3 's. This is because each term contributing to H d in equation (3.5) contains exactly one off-diagonal σ matrix in the tensor product. Therefore, which proves the sublattice symmetry of the Hamiltonian.

Magnetic inversion symmetry
Since the field strength is constant in space and is a two-form, it is invariant under inversion 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 where S ρ is defined as in equation ( Now equations (4.3) and (4.5) fix A d unambiguously. We can write A d explicitly as a matrix through the recursion relation: Note that A d is a Hermitian anti-diagonal matrix. By induction we easily check that A 2 d = 1 so that its eigenvalues can only be ±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 (4.3) as a translation mod 2 along the longest body diagonal of the hypercube: 1 − x k = (1 + x k ) mod 2, for all k = 1, 2, . . . , d. (4.8) Then the position-dependent phase factor in equation (4.5) is exactly the standard phase factor for the corresponding magnetic translation [47].
We now prove A d is indeed a symmetry by induction. The commutator [A d , H d ] is given by By induction hypothesis [A d−1 , H d−1 ] = 0 which is satisfied for d = 2 because H 1 = A 1 = σ 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: where for the last equality we used equation (2.9) and this completes the proof.
Since the symmetry operator is an anti-diagonal matrix, an orthogonal set of eigenvectors is given by (0, · · · , 0, b k , 0, . . . , 0, ±b 2 d +1−k , 0, · · · , 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 ±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 × 2 matrix then it is clear from equation (4.7) that the magnetic inversion can be written as We can check that  (3.5), and the second term becomes but we can freely re-index the summation as because both x k and x c k 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 µν , 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 .

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 (2.3) has only nearest neighbor hoppings. The sum rules can be evaluated in the tensor product representation (3.5) as well. Since this calculation for TrH 4 and TrH 6 in tensor product representation is rather lengthy, we have moved it to appendix A.

TrAH 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 TrH p , what we really should be interested in is 1 2 Tr ((1 + A)H p ). However in this section we will see So for low moments we might as well just study TrH 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 (4.11)), while each factor contributing the Hamiltonian (3.5) contains only one off-diagonal Pauli matrix. It is clear that TrA d H p d = 0 for p < d, because for p < d every term in H p will have at least one diagonal 2 × 2 matrix in the tensor product.

TrH 2
We wish to prove In fact, we will prove a stronger identity for the diagonal entries of H 2 :

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 (5.3).

Tensor product picture
The tensor products in the Hamiltonian (3.5) involve both diagonal and off-diagonal two by two matrices. To contribute to (H 2 d ) ii , terms with the off-diagonal Pauli matrix must be in the same position in both factors of H, so and inside the second term, we have terms It is not hard to see then the sum over x 1 , . . . , x d−1 results in a tensor product of d identity matrices. The same argument applies to (σ 0 ⊗ H d−1 ) 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 Figure 1. Two lattice paths that connect sites i and j with a π flux.

H 2 at π flux
We just demonstrated that the diagonal entries of (H 2 d ) ii = d. We shall further show that at π flux, Together with the sublattice symmetry described in section 4.1, this implies at π flux H d has exactly half of its eigenvalues being − √ d and and the other half being √ d.

Geometric picture
We only need to show Note H 2 d ij receives contributions from 2-step lattice paths that connect lattice sites i and j. 2 There are two scenarios for i = j: 1. There is no 2-step path from i to j. For such pairs of ij clearly (H 2 d ) ij = 0.
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 π, and this means the two original paths give contributions that differ by a factor of e iπ = −1, so their sum vanishes.

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 φ = π. One such pair of cross terms is an anticommutator where the first factor is part of σ 0 ⊗H d−1 and the second factor is the second term in equation (3.5). Since ρ 0 ρ 1 = 0, the product is only nonzero when So the sum reduces to 3 where we have used Note only the first tensor factor in (5.12) depends on x 1 . Now we see the sum over x 1 already gives zero because The same argument can be applied to all other mixed terms. For φ = π we thus demonstrated H 2 d = d1.

TrH 4
We wish to prove TrH 4 does not depend on disorder realizations of S µν .

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φSµν . However, since each clockwise loop is accompanied by its counterclockwise counterpart, the contributions must be functions of cos (φS µν ) = cos φ. 4 We see in both cases the contributions do not depend on the disorder realization of S µν .

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 [27].

TrH 6
In this section we prove that TrH 6 does not depend on disorder realizations of S µν .

Geometric picture
A six-step loop can at most traverse three different dimensions. Let us first think about d = 3. As Parisi argued [27], in three dimensions the field strength tensor φS µν 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 µν . The loops that contribute to TrH 6 can traverse one, two or three different dimensions. Those that traverse one and two dimensions are independent of S µν 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 µν . 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 3cubes they reside in. By the argument just laid out, the sum of each group of such loops is independent of S µν , 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 TrH 6 arises because we have just enough terms in the expansion of TrH 6 to cancel the sine-dependent terms of the form sin φS µν . 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 TrH 8 .

Chaos on the hypercube
For φ = 0 the model is integrable, and has a degenerate spectrum (3.10). The degeneracies are lifted at nonzero φ, but the eigenvalues will eventually flow to ± √ d at φ = π, as predicted by equation ( for the Maldacena-Qi model at infinite coupling. The degeneracies are lifted at nonzero φ, and at φ = π/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 φ < π/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 φ = π/4) which will be studied in more detail below. The apparent crossings of the spectral flow lines are actually avoided crossings even though some are extremely close.

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 [15,40,42,43,48], and can be summarized as follows (see appendix B for more details). Traces of the lattice Hamiltonian, 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 φ. 5 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: where a k is the number of chord diagrams with k crossings. We have defined M HC 2p as a reduced moment since we used TrH 2 in the denominator, but we will call M HC 2p "moment" when the context is free of confusion. In appendix B we lay out the arguments and derivations that lead to equation (6.2) in more details, and discuss the subleading corrections.
The moments given in equation (6.2) are the moments of the density function of the Q-Hermite polynomials: with e 0 = 2/ √ 1 − Q and Q = q = cos φ. However, to include some of the finite-d corrections we set Q = η, which is a renormalized version of q, obtained by matching the fourth moment of ρ QH (x) and the fourth moment of the hypercube model exactly: It is clear η → q = cos φ in the large d limit. In figure 3 we show the average spectral densities for three different values of φ and compare the result with the Q-Hermite spectral density with Q = η. Renormalizing q to η 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 [28] here for the exact results up to the eighth moment: In terms of ρ QH η (x) the spectral density (before ensemble averaging) can be expanded as  Figure 3. The spectral density of the hyper cubic model (black curve) compared to the eighth order Q-Hermite spectral density defined in (6.8) for three different flux values as given in the legend of the figures. In the left figure φ = 0.57π/2, the curve resulted from a fitted a 8 (red curve) differ from the one from the a 8 that is calculated by equation (6.13) (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 φ = 0.95π/2. The deviations are larger for φ = 3 4 π in the right figure.
nomials defined by the recursion relation [49] H with the initial conditions The Q-Hermite polynomials satisfy the orthogonality relation where n η ! is the Q-factorial defined as Note that for the choice of η in (6.4) the coefficients of H η 2 (x) and H η 4 (x) vanish since ρ QH η already gives the exact results for M HC 2 and M 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 TrH 2 and TrH 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) ), where we note again a 6 = a 6 because TrH 6 is independent of disorder realizations, which is not true for a 8 . This is not a good expansion for negative η when n η ! becomes small, see table 1. For example, the large value of a 8 for φ = 3π/4 is due to the smallness of 8 η ! ≈ 0.01. The expansion diverges for φ → π. The reason is that so that a 8 diverges as 1/(q + 1) 2 for q → −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 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 (6.13) but with M HC 8 replaced by the eighth moment of that realization.

Spectral correlations
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 [17][18][19][20][21] 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 π 4 < φ < π/2. Indeed for φ = 0.57π/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 Σ 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 ≈ 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   containing n levels on average is thus δn/n ∼ O(1/d) resulting in a correction to the number variance that behaves as n 2 /d 2 . The results for φ = 0.57π/2 are significantly closer to the random matrix result than those for the other values of φ. For the first (φ = 0.57π/2) and second row (φ = 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 (φ = 3 4 π) a third order polynomial fit to the ensemble average of the spectral density was used to unfold the bulk of the spectrum.  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 difference between the results for ensemble unfolding and local unfolding is due to the fluctuations of a 8 = a 8 + δa 8 . Table 1 contains the results for the simulation parameters of the above figures. We conclude that for φ = 0.57π/2 the collective fluctuations only contribute negligible to the spectral fluctuations, while they are important for φ = arccos(1/13) and φ = 3π/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 φ = 0.57π/2 and φ = arccos(1/13); for φ = 3π/4, where the range of the spectrum where the spectrum can be reliably unfolded is smaller, the width is taken to be 500. For φ = 0.57π/2 local unfolding and ensemble unfolding give almost identical results (see upper right figure of figure 6), while for the other values of φ 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 φ = 0.57π/2 than larger values of φ. 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 [18], 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 [50] 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  Figure 6. The ensemble-unfolded spectral form factor for d = 14 at φ = 0.57π/2, φ = 0.95π/2 and φ = 3π/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 τ 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 φ = 0.57π/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. density is defined by [52]  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 [53]. In the case of hypercube model N = 2 d and ρ(λ) = ρ HC (λ) . In figure 7 we show the microscopic spectral density for an ensemble of 10,000 Hamiltonians for d = 12 and φ = 0.57π/2 (black dots). The result is compared with the analytical result for the chGUE microscopic spectral density (red curve) [54]: where J n (E) are the Bessel functions. We remark that there is no fitting and the agreement is excellent. 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 ±λ k are also visible in the short time behavior of the form factor. Instead of K c (t) ∼ t 2 for the GUE we have K c (t) ∼ t 4 for the chGUE, both for t → 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 Σ 2 (n) ∼ log n) for intervals that are symmetric about zero [55]. 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.

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 [56], and applications of the TFD state can be found in [57,58].
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 (3.8). We thus have where the gamma matricesγ L(R) k are in a representation that was used in [46] to prove that the ground state of the Maldacena-Qi model is a TFD state. Specifically, γ L k =γ k ⊗ 1, k = 1, 2, . . . , d/2, γ R k =γ c ⊗γ k , k = 1, 2, . . . , d/2, whereγ k are Dirac matrices in d/2 dimensions andγ 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γ L andγ R matrices can be obtained by a permutation of the γ L and γ R matrices in equation (3.8) as follows:γ for k = 1, 2, . . . , d/4. Then we can check theγ k matrices in equation (7.3) take the form: for k = 1, 2, . . . , d/4, andγ Since both {γ L k , γ R k } and {γ L k ,γ R k } are Hermitian representations of the Clifford algebra in even dimensions, the similarity transformation U in equation (7.1) that relates the two is unitary. In the Maldacena-Qi model, the basis of the TFD state is constructed from the Hamiltonian is a product of q differentγ L(R) matrices, α is the set of q indices of these gamma matrices, and J α is the Gaussian-random coupling. 6 It is important that the left and right Hamiltonian share the same coupling J α . Because of the tensor structure of the Hamiltonian it is clear that the eigenstates of this Hamiltonian are given by |m ⊗ |n (7.8) with eigenvalues E m + E n . Here, |m are the eigenstates of H L projected onto the the left space. In this basis, the thermofield double state at inverse temperature β is given by iγc Km (7.9) with C R the charge conjugation matrix, and K the complex conjugation operator. In a convention where gamma matricesγ 2k are purely imaginary while theγ 2k−1 are purely real like in equation (7.5), we have that The argument to show that the ground state of the Hamiltonian H M Q is given by the TFD state at β = 0 does not depend on the details of the Hamiltonian (7.7) that determines the basis states [46], 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 mn m|γ kγc |n n|γ kγc |m , (7.12) where in going from the second line to the third line, we have used the fact that 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 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 (7.15) to determine the basis states entering the the TDF state and consider the overlap with the ground state of The gamma matrices in both Hamiltonians are in the representation (3.8). Since the overlap between states is invariant under a unitary transformation, we can do the unitary transformation U in equation (7.1) to transform the Hamiltonians (7.15) and (7.16) into the Hamiltonians (7.7) and the coupling matrix in the right-hand side of (7.2), respectively. Using the above argument, the ground state of (7.16) is given by with eigenvalue E m , then C R K|m is an eigenstate of H R SY K with eigenvalue −E m . The ground state of (7.16) is thus a linear combination of the zero energy states of (7.15). In figure 8 we show the magnitude of the overlap of the ground state with the |m |C R e π 4 iγc Km (denoted by |m | − m 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.
There are other possibilities to choose a basis for a TFD state. For example at zero flux, the Hamiltonian may be written as (7.18) and a TFD state can be constructed out of the eigenstates of H d/2 . For φ = 0 the Hamiltonian has its nonzero matrix elements in the same position as the ones of H d (φ), and also at φ = 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.

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 (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 ±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 δn/n ∼ 1/d, and hence the number variance for large n behaves as Σ 2 (n) ∼ 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 ∼ 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 ±λ 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-toone 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.

A Disorder independence of the fourth and sixth moments in the tensor product representation
A.1 TrH 4 In this section we calculate TrH 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, σ 0 is the 2 × 2 identity matrix, and we refer to equation (3.2) for the definition ρ k .) so that the Hamiltonian of the hypercubic model is given by (see equation (3.5)) The fourth moment can be expressed as Since each h d,µ has only one off-diagonal 2 × 2 matrix in the tensor product, and its position is labeled by µ, the only nonzero traces are of the forms  and in the second line we have used that y k = x − µ + k for k = 1, · · · , ν − 2 so that where we have used σ µ, y( x) ρ x µ−ν = ρ x c µ−ν σ µ, y( x) . It is clear the only nonzero terms are those with Under this condition we see (note that σ ν, y( x) does not depend on the last index of y( x)) whose trace is Trσ µ, x σ µ, x = 2 cos(φS µ,µ−ν ) = 2 cos φ. which is independent of S µν . Combining equations (A.4), (A.5) and (A.14), we obtain the total fourth moment 15) and the normalized fourth moment is equal to which is in agreement with the averaged fourth moment TrH 4 first obtained in [28].

A.2 TrH 6
In this section we show that TrH 6 does not depend on the disorder realization. Most contribution to TrH 6 can be reduced to combinations occurring in TrH 4 . We have two new combinations: Trh µ h ν h µ h ω h ν h ω and Trh µ h ν h ω h µ h ν h ω with µ > ν > ω. For notational clarity we focus on the cases of Trh so that we only need to use with the · · · representing d − 3 additional factors ρ x . Each h appears two times in the traces, and we use another set of dummy indices y, y and y 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 . . 20) and the summation symbol simplifies accordingly: We now work out . The nontrivial part of the trace is where in the last equality we used σ d−2, x σ d−2, y = σ Tr σ d, x σ d, y =2 cos(φS d,1 ) Taking the trace over the remaining ρ x 3 ⊗ · · · ρ x d−1 and sum over the remaining indices x 2 , x 3 , . . . , x d−1 , we finally arrive at which is independent of disorder realizations.
Next we proceed to The nontrivial part is given by Note only the last factor is different from that of equation (A. 22), and the extra conditions enforced this time are So equation (A.26) reduces to Taking the sum of (A.28) and (A.29) we notice the sine terms cancel, thus the result no longer depends on S µν . After performing the sum over x 3 , . . . , x d−1 , we conclude which readily generalizes to the generic cases

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 [27] 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 TrH 2p is given by the sum of all 2p-step Wilson loops: 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) · · · (d − k + 1) ∼ d k . Since k ≤ p, the leading large d contributions will come from those loops that traverse p different dimensions, having a multiplicity factor of d(d − 1) · · · (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. 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 d 2 coordinate axis planes, and if the projection into the µν plane has a plaquette shape, we pick up phase of e iφSµν ; 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φSµν from the µν plane, and the disorder average over S µν on each face results in a cos φ for each projection that loops around a face, and contributions from different plaquettes multiply. Hence we have the following formula: with q := cos φ, where A µν (C) is the area of the loop's projection into µν 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 µ and ν dimensions are represented by the letters a and b, respectively. To study the projection into µν plane, we can temporarily forget letters other than a and b. With regard to a and b, there are only three scenarios: It is clear that the first two cases have zero-area projections in the µν 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 (B.1) and (B.2) for the leading contributions as where q = cos φ 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 (B.7) also appear in the Sachdev-Ye-Kitaev model [15,40,48,59]. The sum on the right-hand side of (B.7) has an interesting solution: in his original paper Parisi [27] 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 [28]. In fact, much earlier this chord diagram sum was studied by Touchard [60] and Riordan [61] in a more combinatorial vein, which led to the Riordan-Touchard formula: 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 (B.7) due to the multiplicity factor d(d − 1) · · · (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) · · · (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 formula 7 it is clear that we can form 1 3 p 2 (2p − 1)!! subleading words of length 2p. 7 We 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 where [m1, m2, . . . , m k ]p denotes a composition of p, that is, an ordered k-tuple (m1, m2, . . . , m k ) such that k l=1 m l = p.
The following map is a bijection between subleading words and interlacing structures of leading words: . where the · · · 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: 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 · · · parts in Eq. B.14, 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. 8 We have demonstrated that the bijection (B.14) 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 (B.14), 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 zeroarea 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: In this case the plaquette projection represented by cdcd cannot be affected because the 8 A 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 1 3 p 2 (2p − 1)!!. We can easily generalize the proof to other intersection structures. Other proofs of this statement already exist, see for example [41,62]. 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  We can obtain a rather compact and visual representation of rule (B.22) 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. represent every chord by a vertex, 2. connect two vertices if and only if the chords they represent intersect each other.
We refer readers to figure 9 for a few examples. In the intersection graph language, the leading moment (B.7) can be written as TrH 2p leading = d(d − 1) · · · (d − p + 1) 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 (B.22), the subleading moment can be written as 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 (B.24) 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 (B.24) 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 leadingcontribution 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 (B.11) and (B.12); the subleading contributions are which are consistent with the results of Marinari, Parisi and Ritort [28]. 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 (B.24), 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 (B.14), namely: . . In some literature [63] 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   otherwise q # → 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.

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. 10 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. . 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.
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: TrH 2p subleading d(d − 1) · · · (d − p + 2) = 1 3 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 G (v 1 ,v 2 ) is the number of edges in G (v 1 ,v 2 ) .