Probing symmetries of quantum many-body systems through gap ratio statistics

The statistics of gap ratios between consecutive energy levels is a widely used tool, in particular in the context of many-body physics, to distinguish between chaotic and integrable systems, described respectively by Gaussian ensembles of random matrices and Poisson statistics. In this work we extend the study of the gap ratio distribution P(r) to the case where discrete symmetries are present. This is important, since in certain situations it may be very impractical, or impossible, to split the model into symmetry sectors, let alone in cases where the symmetry is not known in the first place. Starting from the known expressions for surmises in the Gaussian ensembles, we derive analytical surmises for random matrices comprised of several independent blocks. We check our formulae against simulations from large random matrices, showing excellent agreement. We then present a large set of applications in many-body physics, ranging from quantum clock models and anyonic chains to periodically-driven spin systems. In all these models the existence of a (sometimes hidden) symmetry can be diagnosed through the study of the spectral gap ratios, and our approach furnishes an efficient way to characterize the number and size of independent symmetry subspaces. We finally discuss the relevance of our analysis for existing results in the literature, as well as its practical usefulness, and point out possible future applications and extensions.


I. INTRODUCTION
Symmetry considerations are an essential part of a physicist's toolbox, with countless applications in all fields of physics, ranging from Noether's theorem, gauge theories or the description of phase transitions [1].Another frequent tool is the use of simplified models, which successfully describe the important features of a physical phenomenon without having to deal with all microscopic details.In this respect, Random Matrix Theory (RMT), which was first initiated to understand the statistical properties of energy levels in complex nuclei [2], is an extremely successful approach which has also impacted various branches in physics [3,4].It then comes as no surprise that symmetry properties are an integral part of RMT: one of the best-known examples is the construction of classical Gaussian ensembles from time-reversal symmetry considerations.Depending on the underlying symmetry of the system considered, it is best described by random matrices belonging to one of the three following ensembles: Gaussian Orthogonal Ensemble (GOE), Gaussian Unitary Ensemble (GUE) and Gaussian Symplectic Ensemble (GSE), whose entries are respectively real, complex or quaternionic random variables.Convenient to the description of Floquet operators are the circular ensembles introduced by Dyson [5]: circular Orthogonal Ensemble (COE), circular Unitary Ensemble (CUE) and circular Symplectic Ensemble (CSE).These ensembles have the same asymptotic level spacing distri-butions as the Gaussian ensembles [4].
The celebrated conjectures of Berry and Tabor [6] and Bohigas, Giannoni and Schmit [7] state that RMT describes the spectral statistics of quantum systems with a chaotic semiclassical limit, whereas Poisson statistics provides a description of systems with a classical integrable limit.These two paradigms serve as reference points to study the transitions between localization and ergodicity, for instance the Anderson transition as a function of disorder [8].Quite crucially, quantum many-body systems, for which there is in general no semiclassical limit, also display the same dichotomy: RMT statistics for chaotic systems and Poisson statistics for quantum integrable systems, including those showing an emergent integrability such as many-body localized systems [9,10].Numerous examples illustrate the usefulness of a RMT analysis of quantum many-body spectra [9,[11][12][13][14].
A universal tool in this respect is the study of the distribution p(s) of level spacings, or gaps, defined as the differences between consecutive energy levels, s i = λ i −λ i−1 , assuming that the mean level density is fixed to unity, i.e. s = 1.RMT offers simple, powerful predictions for the distribution p(s) in terms of three different Wigner surmises corresponding to the three Gaussian ensembles mentioned above [2].These surmises are obtained by a simple calculation on random 2 × 2 matrices, and turn out to reproduce most of the features of much larger random matrices, with high precision [15].Normalizing the level spacing distribution requires the knowledge of the density of states, which is often not analytically available.Numerically, one needs to perform an unfolding of the spectrum, for which there exists different procedures [6,13,16].Unfolding can lead to spurious results [17], in particular because of finite-size effects; one may even find instances where different unfolding pro-arXiv:2008.11173v2[cond-mat.dis-nn]25 Feb 2022 cedures leads to different physical interpretations of the same data.For many-body systems, the density of states is generically far from being uniform, which makes the use of the unfolding procedure rather inaccurate.
A very useful alternative to the study of s has been proposed by Oganesyan and Huse [14], in terms of the gap ratio for three consecutive levels, r i = min (si,si+1)  max(si,si+1) .The key point is that considering the ratio of gaps rather than the gaps themselves suppresses the need to know or estimate the density of states, and thus avoids the numerical unfolding step.The probability distribution P (r) of gap ratios is thus well-suited to characterize statistical properties of many-body spectra.The Poisson statistics distribution P Poisson (r) = 2/(1 + r) 2 can be easily derived from a Poisson sequence.For the random matrix spectra, analytical surmises of P GOE (r), P GUE (r) and P GSE (r) have been obtained by Atas et al. [18] from the joint eigenvalue distribution of 3 × 3 random matrices, and improved estimates were obtained in [19] based on 4 × 4 matrices.
Because of its computational advantage (no unfolding needed) and the existence of these analytical predictions, the gap ratio r, in particular its average r and its distribution P (r), has become one of the most studied metrics in the field of disordered quantum systems.For instance, it is often used to characterize the change of statistics across a many-body localization (MBL) transition, between an ergodic phase, for which the RMT predictions for P (r) are expected, and a Many-Body Localized phase, which displays emergent integrability and thus P Poisson (r) gap ratio statistics [9,20,21].The agreement between the RMT-predicted P (r) and the numerical estimate for a given model now routinely diagnoses quantum chaotic models.Any discrepancy in the gap ratio as a function of a model parameter is often interpreted as a sign of a different physical behavior (see e.g.[22]).The distribution of gap ratios is also instrumental in analyzing the symmetry properties of the SYK model and variants as a function of the number of Majorana fermions [23][24][25][26][27]. P (r) has also been measured experimentally to probe an ergodic to MBL transition / crossover [28].Applications of this metrics were also performed in other fields of study, such as in astrophysics [29], for statistics of the zeros of the Riemann zeta function [18], or characterizing entanglement in quantum circuits [30].The computation of the gap ratio statistics has been extended in several ways, such as ratios of gaps for levels with one or more other levels in-between, or non-Hermitian matrices [19,[31][32][33][34][35][36][37].
What happens to spectral statistics in the situation where symmetries are present in the original Hamiltonian?In a seminal work [38], Rosenzweig and Porter computed the level spacing distribution P (s) of systems with several independent random blocks (each being a random matrix with spacing distribution p(s)).This situation typically occurs when a physical system displays discrete symmetries, in which case the number of blocks remains finite in the thermodynamic limit.For contin-uous symmetries, the number of blocks grows with the system size, and ultimately, as many independent spectra are mixed, one expects a Poisson distribution to emerge for large enough systems.In an extension of the original work [38], Berry and Robnik considered mixed phase spaces with both ergodic and integrable blocks [39].
In general, one would be inclined to resolve the underlying symmetries by treating each block independently and performing a block diagonalization.This is not always possible.First, there are cases where a symmetry not previously known or analyzed is discovered fortuitously (e.g. by monitoring the gap ratio and seeing that it does not converge to its expected value).Second, in some situations, the block diagonalization is not known, too complex to implement, or cannot be performed totally.The latter case occurs for instance in systems with non-Abelian discrete symmetries, where two symmetry operations that commute with the Hamiltonian do not commute with each other (see e.g.[40,41]).Third, there are cases where the basis transformation leading to a block structure in the Hamiltonian is known, but results in a Hamiltonian which is more costly to analyze; this is for instance the case if a sparse Hamiltonian leads to non-sparse blocks, inducing a strong decrease in the performances of numerical routines (we will present such an example in Sec.IV D).
In this work, we extend the Rosenzweig-Porter analysis to the computation of the gap ratio statistics P (r) when several independent blocks are present.We do so by calculating the joint gap distribution P (x, y) for a matrix with several independent random blocks, each being a random matrix with joint gap distribution p(s, t).We obtain closed expressions for P (x, y) and P (r) in terms of p(s, t) and its primitives.These expressions are valid for an arbitrary number of blocks.In the case of Gaussian random matrices, we use for p(s, t) a surmise given by the exact 3 × 3 distribution of RMT, which allows us to obtain expressions for P (r).However our formula applies for an arbitrary distribution p(s, t).We note that a recent work [42] provides estimates for P (r) and r based on an surmise obtained from explicit analytical calculations for small-size matrices; however it does not take into account all possible level partitions.Our approach is quite different, as we discuss below.
Our analytical estimates are virtually indistinguishable from numerical simulations on large random matrices.Our results explain several deviations for the distribution P (r) or expectation value r observed in the literature, as discussed in Sec.V A. They can also be useful in several situations such as those mentioned above (which we illustrate with various applications taken from many-body physics in Sec.IV), as well as to estimate the number of effective ergodic blocks in an incompletely thermalized system.
The manuscript is organized as follows.We first introduce the problem in Sec.II, setting up the notations and summarizing the useful literature as well as our own results.Sec.III contains the derivation of the generic form of P (r) when several independent blocks are present.We then present results for the three Gaussian ensembles.Sec.III D compares these analytically obtained results to simulations performed on random matrices, showing an excellent agreement.Sec.IV contains several realistic applications of these results in many-body physics, with a panel of different types of possible symmetries: clock symmetries, symmetries in disorder realizations, dynamical symmetries in Floquet systems, disordered anyonic chains with topological symmetries.In Sec.V, we finally conclude by first discussing existing examples where our work directly applies, and then suggesting some further perspectives.

A. Random matrix ensembles
Let us first consider the case of a single Gaussian random matrix H of size N , whose distribution is proportional to exp(− 1 2 tr H 2 ).We denote by λ 1 ≤ . . .≤ λ N the eigenvalues of such a matrix.The density of eigenvalues is given by the Wigner semicircle law ρ(λ) = [43].Since there are N ρ(λ)δλ levels in an interval δλ, the corresponding mean level spacing in the vicinity of λ = 0 is ∆ = π/ √ 2N , which gives a local density 1/∆ = √ 2N /π.The joint distribution of eigenvalues is [4] β = 1, 2 or 4 is the Dyson index and N and a are normalization constants.
In a region of constant density, the nearest-neighbour spacing distribution p(s) is well-approximated by the Wigner surmise [2], corresponding to the exact result obtained from Eq. (1) for 2 × 2 matrices, where a β , b β are normalization constants, chosen in such a way that s = 1.In a similar way [18], one can approximate the joint distribution of consecutive nearestneighbour spacings p(s, t) by its exact expression for 3×3 matrices, which can be obtained from Eq. ( 1) with N = 3 by integrating over one variable.It reads where the constant B β is such that both spacings are normalized as s = t = 1, and A β is the overall normalization factor.From this expression, one can then obtain the distribution of r = min(t/s, s/t) as Since the distribution p(s, t) is symmetric in s and t, Eq. ( 4) reduces to This approach was carried out in [18], yielding with Z β the normalization constant.Since Gaussian and circular ensembles have the same asymptotic level spacing distribution, the same analysis should equally be valid for circular ensembles, the only difference being that the mean level spacing is ∆ = 2π/N and thus the density of states is uniform and proportional to N , rather than circular and proportional to √ N .For finite N , this difference in the shape of the density can result in small differences between the circular and Gaussian ensembles which are expected to vanish in the largematrix limit.For the 3×3 matrices leading to the surmise Eq. ( 7), the difference is already very small [44].

B. Compound spectrum: the Rosenzweig-Porter approach
Let us now consider ensembles of random matrices of size N which can be decomposed into m independent blocks of sizes such a matrix can be obtained by diagonalizing each block separately and ordering the eigenvalues, so that the spectra of the blocks are interlaced.Let N = (N 1 , N 2 . . .N m ) be a vector of block sizes.The compound spectrum {λ i , 1 ≤ i ≤ N } can be characterized by its spacing distribution P N (s), which is the distribution of gaps s i = λ i+1 − λ i .It can also be characterized by the gap ratio distribution P N (r), with ri = s i /s i−1 , or [45] by the gap ratio distribution P N (r), with If there is a statistical symmetry between left and right intervals then the relation P N (r) = 1 r2 P N ( 1 r ) holds, which entails that P N (r) = 2P N (r)θ(1 − r) [18].In that case, the distributions of r and r essentially contain the same information.As we shall see, this is the case for the distributions considered in this paper, and therefore, as is often done in numerical simulations, we concentrate on the distribution P N (r) with r ∈ [0, 1].
If the m blocks are independent Gaussian random matrices given by the Wigner-Dyson ensembles with index β, then the spectrum of block i, {λ characterized by its mean level spacing around λ = 0, given by ∆ i = π/ √ 2N i , or by its local density ρ i = √ 2N i /π.The resulting spectrum obtained by the superposition of the m spectra has density ρ = i ρ i .Introducing the normalized densities µ i = ρ i /ρ, we have µ i = N i /N .For the circular ensembles, where densities are uniform over the unit circle, µ i = N i /N .The Rosenzweig-Porter approach, which gives the nearest-neighbour spacing distribution P (x) associated with the compound spectrum, consists in assuming that the compound spectrum is a superposition of independent and identically distributed spectra with uniform density ρ i and with nearest-neighbour spacing distribution given by the surmise Eq. ( 2).The computation proceeds by identifying that a gap in the compound spectrum can originate either from a gap in one of the spectra or from a gap between eigenvalues from two distinct spectra.Considering all possibilities and the probabilities attached to them leads to the spacing distribution P (x).This approach is detailed in Sec.III A. These results were extended in [39] by considering a mixed phase space, which amounts to adding Poisson blocks to a chaotic Hamiltonian.The work [39] derives an explicit formula for P (x) for a Poisson block and (N − 1) chaotic blocks with same density.

C. Summary of our results
In this paper we extend the Rosenzweig-Porter approach to derive the joint distribution of consecutive nearest-neighbour spacings P (x, y) of a compound spectrum made out of several spectra with arbitrary distribution.It is given by the very compact expression Eq. ( 27)- (28), for which we give a probabilistic interpretation.We then obtain P (r) from the analog of Eq. ( 6), namely Applying our expressions to the RMT expressions Eqs.
(2) and (3), we obtain a closed general expression for P N (r).We then apply this general formula to the case of identical block sizes N i = N/m, for which we use the short notation P m (r).Some of these calculations result in exact closed (albeit complex) forms, others require a numerical integration.Besides the full distribution, we will also consider the average gap ratio r m = 1 0 rP m (r)dr and the limiting value for vanishing gap ratio P m (0) = lim r→0 P m (r), as they turn out to be of great practical use to identify the existence of a symmetry (P (0) = 0 in the no-symmetry case m = 1).In Section III E we also consider the quantity I 1/4 m = 1/4 0 dr P (r), which proves useful to identify symmetries in an experimental setting where few realizations of the spectrum are available.Our results are summarised in Table I.In the electronic supplementary material, we provide a Mathematica notebook allowing to reproduce our calculations.

III. ANALYTICAL RESULTS
We now turn to the detailed proofs of our analytical formulae.The reader not interested in the details of the Since the function p(s) in Eq. ( 2) corresponds to spectra with mean level spacing equal to 1, the spacing distributions of each spectrum are given by the function Eq. ( 2) rescaled by the mean level spacing, i.e. p(ρ i s).We introduce the functions and The function f (s) gives the probability to have λ i+1 ≥ s knowing that λ i = 0.The function g(s) gives the probability to have λ i+1 ≥ s knowing that λ i ≤ 0, that is, the probability to have a spacing at least s.These probabilities are related through the identities Introducing the rescaled spacing x = ρs we have from Eq. ( 11) Spacings arise as empty intervals ]λ q [ of length s with i, j = 1, ..., m.We have to consider the two possibilities i = j or i = j, and calculate the probability densities associated with each configuration.
Let us first consider the case where m = 2 (we will later generalize this analysis to more blocks).We have to consider the two following cases: 1. Configurations giving rise to an empty interval of type ]λ q+1 [, which are are such that λ q +1 for some q .The probability of such a configuration for i is given by p(ρ i s), while the probability for j is g(ρ j s) since λ (j) q and λ (j) q +1 can be anywhere outside ]λ Taking into account the probability µ 2 i to have a level i at both ends of the interval, we get for the configuration ]λ q+1 [ a probability density µ 2 i p(ρ i s)g(ρ j s).Using Eq. ( 12) we can rewrite it as [∂ 2 x g(µ i x)]g(µ j x). 2. Configurations giving rise to an empty interval of type ]λ q [, which are such that λ q+1 .The probability of such a spacing for i is given by f (ρ i s) since λ (i) q+1 can be anywhere in ]λ (j) q , ∞[, while the probability for j is f (ρ j s) since λ (j) q −1 can be anywhere in ] − ∞, λ (i) q [.The probability of having a level i and a level j at the ends of the interval is given by µ i µ j , so that the probability density of configuration ]λ Summing these probabilities over i, j = 1, 2 we get for the spacing probability P (x) of the mixed levels or equivalently The above reasoning proves Eq. ( 14) for m = 2.In fact, given the final expression, we can come up with a much shorter proof of the validity of Eq. ( 14) for arbitrary m.Indeed, the probability of finding an interval of a given length between two consecutive eigenvalues is the second derivative of the probability of finding an interval larger or equal to it with no eigenvalue in it (see Eqs. (10) and ( 11)).Therefore P (x) must be the second derivative of the probability of finding an empty interval larger or equal to x.The probability that no level of type i occurs in an interval of size x is g(µ i x).Since levels from different sequences are independent, the probability that no level of any type occur in an interval of size x is simply the product of all g(µ i x), which directly entails Eq. (14).Incidentally one can check, using Eqs.( 9)- (11), that P (x) in Eq. ( 14) is properly normalized to 1.

B. Joint consecutive spacing distribution P (x, y)
We now apply the same line of reasoning to the joint distribution of two consecutive spacings.Our aim is to obtain the joint distribution P (x, y) in terms of the distribution p(s, t) for a single spectrum.
Starting with the function p(s, t) in Eq. ( 3), we introduce the function which is the marginal distribution of p(s, t).Since p(s, t) defined in Eq. ( 3) is symmetric in the exchange of s and t, the marginal distribution can be equivalently taken by integrating over the first variable.Note that this expression differs from the one in Eq. ( 2), which corresponds to the result for 2 × 2 matrices while p(s, t) was obtained for 3 × 3 matrices.Functions f and g can then be defined from p as in Eq. ( 9)- (10).In terms of p(s, t), their explicit form is and We also define the two-variable functions These functions are related through the identities The four configurations of consecutive spacings considered in Sec.III B. Different colours correspond to distinct spectra.The three central levels are the ones from which the ratios are calculated; the outer levels are at the same height to indicate that their relative position is irrelevant.
The analogs of Eqs.(12) are A sequence of two consecutive spacings arises as a sequence of two empty intervals ]λ q , [ with i, j, k = 1, ..., m.We have to consider all possibilities for i, j, k and calculate the probability densities associated with each configuration.
Once again, let us consider the simplest case m = 3, that we will later generalize.We only need to examine four cases, corresponding to patterns iii, iij, iji and ijk and depicted in Fig. 1 : q+2 , which arise whenever λ q +1 for some q and λ (k) q +1 for some q .The probability of such a configuration for i is given by p(ρ i s, ρ i t); the probability for j is g(ρ j (s + t)) since λ (j) q and λ (j) q +1 can be anywhere outside ]λ q+2 [; and the same goes for k.Taking into account the probability µ 3 i to have three levels i at the ends of the intervals, we get for this configuration a probability density µ 3 i p(ρ i s, ρ i t)g(ρ j (s + t))g(ρ k (s+t)).Using Eqs. ( 21)-( 24) we can rewrite it as

Configurations λ
q , which arise whenever λ q +1 for some q .The probability for i is e 2 (ρ i s, ρ i t), the probability for j is f (ρ j (s + t)), while the probability for k is g(ρ k (s + t)).We get for this configuration a probability density We used the fact that ∂ y g(x + y) = g (x + y).

Configurations λ
q +1 for some q .The probability for i is p(ρ i (s + t)), the probability for j is h(ρ j s, ρ j t), and the probability for k is g(ρ k (s + t)).We get for this configuration a probability density We used the fact that ∂ x ∂ y g(x + y) = g (x + y).

Finally, configurations λ
q , which arise whenever λ and λ q+1 .The probability for i is f (ρ i (s + t)), the probability for j is h(ρ j s, ρ j t), and the probability for k is f (ρ k (s+t)).We get for this configuration a probability density We can now sum all contributions over i, j, k = 1, 2, 3.There are 27 terms, which can be put under the compact form Equation ( 25) has a simple probabilistic interpretation.
Let us define a function H as by analogy with Eq. (19).Thus H(x, y) gives the probability of having a triplet (λ q−1 , λ q , λ q+1 ) of levels of the mixed spectrum such that λ q−1 < λ q − x and λ q + y < λ q+1 .That is, H(x, y) is the probability that some level λ is such that all other levels λ (j) q verify either λ (including the case i = j, in which case of course q = q).At fixed i, the probability of such a configuration is h(µ i x, µ i y) for spectrum i, and g(µ j (x + y)) for all j = i.Summing over all i (and taking into account the probability µ i to have a level i in the middle), we get for H(x, y) the expression under the derivation symbols in Eq. (25).
In fact, this reasoning provides a proof of the general case with arbitrary number m of spectra.We thus have in the general case with One can check, using Eqs.( 20)-( 24), that P (x, y) in Eq. ( 27) is properly normalized to 1.Note that, although in what follows we will apply Eqs. ( 27)-( 28) to the random matrix case, these equations are valid for an arbitrary initial distribution p(s, t) of individual spectra.
C. Ratio distribution P (r) In order to obtain P (r), one first needs to evaluate functions g and h to obtain H(x, y) using Eq. ( 28), then take its derivative with respect to x and y, and finally perform the integration in Eq. ( 8).In the GUE and GSE case, there is no closed-form expressions for the function h, so that we are left with a double integral (one in the definition of h, one corresponding to the final integration in Eq. ( 8)).In the GOE case however, we obtain explicit expressions for g and h, as we will show below, and thus we get a closed-form expression for P (x, y).The remaining integral Eq. ( 8) is doable analytically only in the case of a mixture of m = 2 spectra.

GOE case
In the GOE case, the joint distribution p(s, t) reads Starting from Eq. ( 29) for p(s, t) and calculating explicitly the functions g and h given in Eqs. ( 17) and ( 19) we get with One can then rewrite Eq. ( 28) as with a1...am some polynomials of x, y and the µ i , which can be obtained explicitly from Eqs. (30) and (32).
Functions U i have the property that they transform into each other under derivation.Namely, the derivative of any function i c i U i (with c i polynomial in s) is of the form i ci U i (with ci polynomial in s).The same property holds for the V i upon derivation with respect to s or t.Therefore, using Eq. ( 27) and the expansion Eq. ( 34), we obtain where a1...am are polynomials of x, y and the µ i .Given the definition of the functions U i and V i , Eq. ( 8) can be expanded as a linear combination (with real coefficients dependent on the µ i ) of integrals of the form with λ and the a i depending on the µ i and on r (and possibly a i = 0).It appears that in general such an integral does not have a closed form.However in the case m = 2 we have the identity from which one can deduce Eq. ( 36) for all odd values of k by deriving with respect to λ, and for all even values of k by first integrating by parts and then deriving with respect to λ.This yields a (rather lengthy) closed-form expression for P (r) in the case m = 2 (which is given in full in the electronic Supplementary Material).To give an idea of this expression, we only give P (0) in the case of two blocks of the same size: For m blocks of the same size, we get P (0) = 1.71587 for m = 3 and P (0) = 1.83279 for m = 4, as reported in Tab.I.In practice, the fastest way of obtaining P (r) for GOE in the general case is to calculate H(x, y) and P (x, y) analytically from the explicit expressions for h and g, using Eq. ( 27)- (28), and perform the last integral numerically.The Mathematica notebook in the electronic Supplementary Material implements the two possibilities to obtain P (r).From the exact equation, numerical integration over [0, 1] yields the mean ratio.For instance for m blocks of equal size, we get r GOE,m blocks = 0.423415, 0.403322, 0.396125 (39) for m = 2, 3 and 4 blocks respectively, again reported in Tab.I.

GUE case
In the GUE case, the joint distribution p(s, t) reads The calculation of g and h for GOE was made possible by the fact that either s or t is of degree 1 in the polynomial in front of the exponential in Eq. ( 29).This is no longer the case for GUE and GSE.However, g as well as the first derivative of h can be obtained analytically.We find where T (x, a) is Owen's T -function, defined as For h we have h(x, y) = ∞ 0 e 2 (x + a, y)da, with and e 1 (x, y) = e 2 (y, x).Using the identities in Eq. ( 20), the expression of P (x, y) reads in which only the first sum involves the function h for which no closed form is available.The ith term in that sum is µ i h(µ i x, µ i y)q i (x + y), with q i an explicitly known function involving only products of derivatives of g.For this term the integral Eq. ( 8) yields a contribution We perform numerically the twofold integrals Eq. ( 45) and the single integrals over x for all the other terms.For m blocks of equal size, we obtain: and r GUE,m blocks = 0.422085, 0.399229, 0.39253 (47) for m = 2, 3, 4 respectively.

GSE case
For GSE the joint distribution p(s, t) reads The function g reads with and It can be checked that the second derivative of g is indeed the marginal probability p, and that g(0) = 1.Similarly as in the GUE case, h is not calculable in closed form, but we have h(x, y) = ∞ 0 e 2 (x + a, y)da with Erf 243 320 where R 1 (x, y) and R 2 (x, y) are polynomials given by and The computation is then the same as for GUE.Using Eqs. ( 44) and ( 45) we get, for m = 2, 3, 4 blocks of equal size, We finally mention that the Mathematica notebook in the electronic Supplementary Material allows to reproduce these computations, as well as to consider different cases (other values of m, unequal sizes of the m blocks).

Poisson (m → ∞) limit
One can easily check that in the case of a mixture of m spectra of the same size, one recovers the Poisson distribution in the m → ∞ limit.Indeed, in that case the function H reads and thus with functions h and g evaluated at x m , rx m and (1+r)x m , respectively.In the limit m → ∞, these arguments go to 0. From the explicit expressions for g and h, the only term in the square brackets that survives is g 2 , which goes to 1. Using the fact that g(x) = 1 − x + O(x 3 ) close to 0, we get which is indeed the Poisson result.

D. Comparison with numerics
We now compare the gap ratio distribution P m (r) and the mean value r m obtained through the analytical approach of Sec.III to direct numerical computations on large random matrices.Numerical RMT spectra are computed using the matrix models of [46], based on tridiagonal matrices.These models are numerically faster to diagonalize, but otherwise equivalent to the dense ones.Full lines are the surmises obtained from Sec. III, except the m = 1 case, for which the corresponding surmise Eq. ( 7) is taken from [18].Points are numerical results from random matrices of size at least Ntot = 2000, averaged over 3.6 × 10 5 histograms.
total linear size.In all numerical computations presented here, N tot is at least 2 × 10 3 , and 3.6 × 10 5 samples are used.
Figure 2 displays the results of this comparison for m = 2, 3, 4 (as well as the surmise of [18] for m = 1) and all Gaussian ensembles.The comparison is excellent and within the scale of this figure, there is no visible difference between the analytically obtained P m (r) and the numerical results P num m (r).More precisely, we found that the relative error |P m (r)−P num m (r)|/P m (r) is always less than 0.01.This translates into also an almost perfect agreement (with no difference within error bars) between r m (from Tab.I) and the numerical estimates reported in Tab.II.
We also compare analytical and numerical results for the specific case of m = 2 GOE blocks of different sizes.In Fig. 3 we present results for the probability distribu-m GOE GUE GSE r [18] 1 0.5307(1) 0.5996(1) 0.6744(1) r m 2 0.4235(5) 0.4220(5) 0.4116(5) 3 0.4035(5) 0.3992(5) 0.3927(5) 4 0.3963(5) 0.3924(5) 0.3886 (5) TABLE II.Values of averages r for m blocks, as obtained from simulations on random matrices.The value for m = 1 is taken from [18].Notice the excellent agreement with the surmise results reported in Tab.I. tion P (r) for different values of α = µ1 µ1+µ2 (top panel) as well as for the expectation value r 2,α as a function of α (bottom panel).Here again the agreement between analytical and numerical results is striking.More precisely, we found that the relative error |P m (r)−P num m (r)|/P m (r) was always less than 0.01 for block ratios α ≥ 0.2.For α < 0.2, the relative error increases, but remains below 0.05.The seemingly crossing point in Fig. 3 is in fact not a crossing point, as one can convince oneself by using the exact expression for P GOE 2 (r) and calculating its value with enough precision in the vicinity of that point.
We conclude this comparison section by discussing [42], which provides an analytical estimate for P (r) for m = 2, derived from a 4 × 4 surmise which is forced to contain two levels of each of the m = 2 blocks.This does not contain all possible patterns considered in Sec.III B. The estimated value for r obtained from this approach is approximately close to the one presented in Table I for the GOE, GUE but strongly differs for the GSE, while our results in this latter case agree with the numerical estimates in Table II.

E. Symmetry detection in an experimental context
The numerical comparisons in the previous subsection are done for matrix sizes N ∼ 2000 and histograms are obtained from many random realizations.This allows to compare to the analytical results, obtained in the thermodynamic limit, with small enough statistical error bars.In the context of numerical simulations, the quantities r m and P m (0) given in Table I provide a signature allowing to identify the presence of a symmetry.This generalizes the existing results for r which is a quantity routinely used in numerical studies (see discussion in Sec.I) to identify the chaotic nature of a spectrum.
While most of the applications of gap ratio statistics have indeed been so far used inputs of numerical spectra of many-body systems, it is worth discussing applications to experimental spectroscopies, which typically involve less statistics (less realizations of disorder) and (sometimes) smaller spectra.For instance, the experimental measurement of P (r) in Ref. [28] were performed in a system with N = 45 energy levels, and using 4 realizations of disorder.Other typical experiments probing disordered many-body quantum systems (often in the context of many-body localization) in various experimental platforms (cold atoms, trapped ions, superconducting qbits) average experimental results over 6 [47][48][49][50][51], 12 [52], 20 [53], 24 [54], 30 [55], 50 [56,57], and up to 197 [58,59] realizations of disorder.Most of these platforms work on quantum systems with a minimum of tens of qbits or atoms with corresponding many-body spectra of at least N = 1000 energy levels.In a different physical context, spectroscopy experiments on nuclei allow to resolve a quite large number of energy levels (often by combining results from different experimental techniques), typically from hundreds to thousands [60,61].
Interestingly, already at sizes achievable experimentally our approach provides a signature of symmetries.As can be seen in Table I, differences between values of r m are quite small.Having experimental investigations in mind, we therefore propose to consider instead the quantity which is simply the integral of the distribution of r up to a point chosen at r = 1/4; this upper bound is arbi- trary, but it is close to the crossing point r ≈ 0.288 of P GOE m=2 (r) and P GOE m=3 (r).One can easily obtain a numerical estimate of I 1/4 m from an experimental spectrum by counting the number of ratios less than 1/4.From the analytical side, theoretical expressions can be obtained from our exact formulas and are given in Table III E.
To illustrate this approach, we give an example of a 'numerical experiment' where one would like to distinguish between the cases m = 2 and m = 3 in a case where the total number N of available levels is small and realizations are scarce.In Fig. 4 we display probability distributions for the quantity I 1/4 m when data are collected from spectra of size N = 180 and when 40 realizations of the experiment are available (solid lines in Fig. 4).In such a case, the number of available levels is very small since each block has size N/m (90 for m = 2 or 60 for m = 3).The two histograms associated with m = 2 and m = 3 are clearly distinguishable.In the case of an even smaller size N = 48, one needs about 120 realizations to get a comparable width of the histograms.These values of N and number of realizations of disorder are comparable to the experimental situations discussed above.
More interesting is the probability of correct identification of the symmetry.If m = 2, the experimentally measured value would be smaller than defined by Eq. ( 60) obtained from the surmise approach in Sec.III.

IV. ILLUSTRATIONS IN QUANTUM MANY-BODY PHYSICS
We now illustrate the usefulness of the above results by comparing them with simulations on realistic spectra obtained from quantum many-body problems.Most of our examples are taken from one-dimensional lattice models, mostly for computational convenience.In the following, the lattice will thus be a one-dimensional chain with L sites.Except otherwise mentioned, we will explicitly break translation symmetry, as well as possibly other lattice symmetries (such as reflection around the center of the chain) to concentrate on the existence of a few blocks.The existence of translation symmetry would result in the existence of L blocks (labeled by the L reciprocal wave-numbers), which would result, as discussed earlier, in an (effective) Poisson distribution for level spacings and gap ratios in the thermodynamic limit.The translation symmetry will be broken by using disorder characterized by a disorder strength .In all the simulations presented below, we take not too small (in order to avoid the proximity to the translation-invariant case, which would cause stronger finite-size effects) as well as not too large, to avoid for instance a possible many-body localized phase (which would also result in Poisson spectral statistics).In all Hamiltonian systems we examine, we consider mid-spectrum eigenstates, obtained either by full diagonalization (for the smaller Hilbert space sizes) or by the shift-invert subset diagonalization method [62] for larger systems.

A. Quantum clock models
The first example deals with Q-states quantum clock models, which are natural Z Q -symmetric generalizations of the Ising quantum chain with Q-states quantum "spins" on each site [63,64].These exhibit a rich ground state phase diagram including ordered and disordered phases as well as critical lines, and have attracted a lot of attention in the recent years due to their relation with parafermions, a Z Q generalization of Majorana fermions [65], as well as with topological phases [66].They are furthermore related to cornerstone models of statistical mechanics, including the Potts model (where the Z Q symmetry is promoted to a larger, S Q symmetry) and the chiral Potts model [67,68].
On each site, we define a spin taking Q possibles values (0 . . .Q − 1), as well as two operators σ and τ , which generalize the Pauli matrices σ z and σ x of the Ising chain: σ measures the orientation of the spin, while τ rotates it by one unit "around the clock", and as a result these fulfill the following algebraic rules: Simple matrix representations are obtained in the basis where σ is diagonal (the "{σ}-basis"): In the basis where τ instead is diagonal (the "{τ }-basis"), the matrices are exchanged.The standard Hamiltonian for quantum clock models is written as a linear combination of (τ j ) a , a = 1, . . ., Q−1 on each site j and exchange terms (σ † j σ j+1 ) a , a = 1, . . ., Q − 1.It is invariant under a Z Q "clock" symmetry σ j → ωσ j , and the associated conserved charge Z = j τ j has eigenvalues 1, ω, . . .ω Q−1 .For Q ≥ 3 the original model has two other important symmetries: charge conjugation, which acts as τ j → τ † j , σ j → σ † j , and time-reversal, which is anti-unitary (and therefore sends any constant to its complex conjugate) and sends σ j to σ † j while leaving τ j invariant.The model we consider in the following breaks all symmetries, but Z Q : where the sum runs over the L sites of the 1d lattice.For practical computations we restrict ourselves to Q = 2, 3, 4. The coupling constants J j are independent random numbers uniformly taken from a box distribution [J − , J + ].Since they are a priori different on each site, they break invariance under translation or spatial reflection.The last term breaks both time-reversal and charge conjugation symmetry (this breaking could also have been achieved by perturbing with the U (1) charge S z introduced in [69]).
We first consider results of simulations performed in the {σ} basis, with a full Hilbert space of size Q L .In the top panel of Fig. 5 we present the average gap ratio for different chain sizes.We clearly observe gap ratios which do not tend to their GOE value r GOE , but FIG. 5. Q-states clock model -Top: Average gap ratio r for the model of Eq. ( 62) for different values of Q, as a function of Hilbert space size |H|.Open symbols denote simulations where the ZQ symmetry is resolved, in which case the Hilbert space size is the block size Q L−1 (we averaged data over all Q equivalent blocks).Filled symbols denote results for the full Hilbert space of size Q L .The dashed lines represent the values of r GOE m obtained in Sec.III (taken from Table II), while r GOE is the numerical estimate for the GOE distribution taken from [18].The precision of our numerics allows us to clearly distinguish the case m = 4 blocks with r GOE m=4 from the Poisson value r Poisson also represented in the plot.Simulations parameters are = 0.5, J = 1, Γ = 0.8, g = 0.5.For Q = 2, instead of the time-reversal and charge conjugation breaking term in Eq. ( 62) which vanishes when Q = 2, we add a next-nearest neighbor interaction g2 j σjσj+2 in order to break the mapping to a free-fermion model (we take g2 = 0.5).Statistics are obtained by focusing on 200Q eigenstates in the middle of the full spectrum, except for the smaller sizes where ∼ 20Q eigenstates where considered.Results are averaged over more than 4000 realizations of disorder, except for the largest size where 1000 realizations were used.For Q = 2, 3, 4, we obtained results on chains of sizes up to L = 17, 11, 9 respectively.Bottom: Probability distribution of the gap ratio P (r), as obtained from simulations of chains of sizes L = 16, 10, 8 for Q = 2, 3, 4 respectively.Simulation parameters are the same as in top panel.The solid lines represent the surmises P GOE m (r) obtained from the analytical computations in Sec.III.rather to their r GOE m=Q value as the size of the Hilbert space is increased.This is expected, as the Hamiltonian possesses Q sectors of identical size Q L−1 labeled with the different eigenvalues of the charge Z.Note however that working in the {σ} basis does not allow to simply construct the Hamiltonian blocks, as Z is off-diagonal in that case.Furthermore, the Hamiltonian is complex in this basis, and without any further indication on the existence of the Z Q symmetry, it would not be clear why GUE statistics should not show up.
When switching to the {τ } basis, Z is now diagonal and the blocks are easily constructed.Furthermore the Hamiltonian becomes real.Computing the average gap ratio in each block leads to an asymptotic r GOE value for each block, showing that each block is indeed independent and no further symmetry has been missed.
We further confirm these results by showing the full distribution P (r) in the bottom panel of Fig. 5 for Q = 2, 3, 4. When the full spectrum is taken, the distribution obtained numerically for the largest system size is in excellent agreement with the surmises P GOE m=Q (r) obtained from Sec. III.
Note the importance of the time-reversal symmetry breaking term g = 0 in Eq. ( 62) in this analysis.In the presence of time-reversal symmetry (at g = 0), the blocks with Z and Z * are identical, leading to exact degeneracies.These additional values at r = 0 would result in effective values of r lower than their Poisson values for finite-size systems.

B. Discrete symmetries in disorder distributions
In this section, we consider the Heisenberg spin chain in the presence of a random external field h j : where σ x,y,z are the standard Pauli matrices, and we use periodic boundary conditions.In general, this system hosts a many-body localized phase at large enough disorder: in particular, the model with box disorder has become the standard model of MBL in one dimension [9,21].When disorder is reduced, the system undergoes a transition towards a thermal phase, a signature of which is an RMT-like spectral statistics.Since uncorrelated disorder explicitly breaks all spatial symmetries, we expect the spectral statistics in the thermal phase to be of GOE type.However, in the specific case of binary disorder, i.e. taking on discrete values h j = ±h, visible deviations from the GOE gap ratio distribution were observed in the bulk of the thermal phase [70].The authors of [70] explained that this phenomenon was due to the peculiarities of discrete disorder distributions.Indeed, while on a finite-size system a typical disorder configuration will break all spatial symmetries, with discrete disorder distributions such as the binary one, there is a non-zero probability that one or several of them are preserved, specially when considering periodic boundary conditions.For example, out of the 2 4 = 16 possible binary disorder configurations on L = 4 sites, 4 are reflection symmetric: (+h, +h, +h, +h), (+h, −h, −h, +h), . . .(actually, all disorder configurations on L = 4 sites have a spatial symmetry: reflection, translation, inversion -exchanging h ↔ −h, or a combination of them).Of course, when L is increased, the probability of drawing a spatially symmetric configuration decreases exponentially fast; but for the largest system sizes within reach of exact diagonalization techniques, the fraction of symmetric disorder configurations is still large enough to significantly alter the level statistics if disorder averaging is done "naively", that is, by uniformly sampling over disorder.A possible workaround put forward in [70] is to discard the spatially symmetric disorder configurations.If one insists on using all samples, another possibility to is to explicitly resolve the symmetry block structure of the Hamiltonian, whenever the disorder configuration happens to be symmetric.This is cumbersome, especially given the number of possible symmetries that must be taken into account.63) with L = 18 spins, with periodic boundary conditions.Data is averaged over all realizations of the binary transverse field hj = ±1/2 (2 18 in total, 3914 of which are nonequivalent up to symmetries), and over 150 eigenstates around infinite temperature energy E = (Emin + Emax)/2.The red solid line represents the analytical prediction Eq. ( 64): a linear combination of surmises P GOE m (r) obtained from the analytical computations in Sec.III.For comparison, the green solid line shows the predicted distribution when two blocks contributions are not taken into account.Inset: Difference δP (r) between the numerical data and these surmises.
In that context, using the surmise for several GOE blocks proves useful.Indeed, the gap ratio in the thermal phase of the model can be written as a sum over symmetry sectors: where w m is the weight of the symmetry sector of m blocks.In the thermodynamic limit L → ∞, w 1 → 1, while w m>1 decays exponentially fast to zero.Note that for samples with two blocks (m = 2) the two blocks are always of the same size, whereas for samples with more than 2 blocks, m > 2, blocks are not necessarily of equal size.While the expressions in the previous section allow us to compute the surmise for these non-homogeneous samples, we can to a good degree of approximation neglect their contribution to the gap ratio distribution.Indeed, we find using simple combinatorics, that the total weight m>2 w m coming from samples with more than two blocks is more than halved when L → L + 2, and for L = 18, it already represents less than 0.2% of the total weight.We will therefore make the approximation that P (r) = w 1 P GOE 1 (r) + w 2 P GOE 2 (r).In Fig. 6, we show the gap ratio distribution for the Hamiltonian Eq. ( 63) for L = 18.We find w 1 = 243936/2 18  0.93, w 2 = 17640/2 18 0.07.This system size is of the order of what is achievable using stateof-the-art exact diagonalization techniques targeting the middle of the energy spectrum [62].However, it is not large enough for w 2 to be negligible compared to w 1 .Indeed, as shown in Fig. 6, incorporating the m = 2 contribution visibly improves the agreement with numerical data.Note the clear difference at r = 0 between Eq. ( 64) (for which P (0) = 0, as in the numerical simulations of Eq. ( 63)) and P GOE (r) which vanishes at r = 0 due to level repulsion.Accordingly, the predicted average gap ratio using the m = 2 surmise r = w 1 r 1 + w 2 r 2 0.527 is closer to the numerically computed value r − r Heisenberg 0.004 than the "naive" prediction involving only m = 1: r m=1 − r Heisenberg 0.013.

C. Floquet spin chain model
We next consider both a static and a Floquet spin 1/2chain model.Floquet systems have attracted a great deal of interest, because while they are amongst the simplest non-equilibrium Hamiltonian systems, they exhibit new non-trivial properties, that are not observed in their static cousins.In particular, single-particle Floquet systems can host topological phases that have no static equivalent [71].Interacting many-body Floquet systems a priori exhibit no such interesting phases of matter, since the combination of interaction and driving is expected to heat up the system to an infinite-temperature, featureless state [72][73][74].However, it has been shown [75][76][77][78] that the addition of disorder, hindering energy propagation throughout the system via a MBL mechanism, can prevent heating and give rise to new interacting Floquet phases, such as discrete time crystals [79,80].Here, we study an interacting Floquet system, along with its static counterpart, for comparison.We show in the following that the Floquet system exhibits an extra symmetry, that can be associated to Floquet topological modes.In order to detect the symmetry, we adjust the system parame-ters so as to be in the thermal phase of the model.Then, level statistics is expected to follow RMT predictions, enabling us to employ our surmises to detect the Floquet symmetry.
We work with the following spin 1/2 Hamiltonians: again with periodic boundary conditions, which we combine to form a time-independent H static = H x + H z and a time-dependent model: Because the drive is periodic H driven (t + τ ) = H driven (t), such a model is indeed a Floquet system.In the Floquet setting, energy is not conserved.It is replaced by quasi-energy, which is defined up to arbitrary shifts by 2π/τ .More specifically, let us introduce the Floquet operator U F = exp(−iτ H x ) exp(−iτ H z ), which is the evolution operator over one drive period.To the unitary Floquet operator we can associate a Floquet Hamiltonian H F , defined as U F = exp(−iτ H F ), whose eigenvalues ε α and associated eigenvectors are respectively the quasi-energies and the Floquet eigenstates, which hold information about the dynamics and steady-state properties of the system [81].In practice, when computing level statistics, we will therefore use the quasi-energies ε α exactly like energies in the static case.
Going back to the system Eq.( 65), remark that the random longitudinal fields h j = 0 break both the Ising and the translation symmetries.Our model differs from the most commonly used one in that h j = 0 on even sites.This does not change the physics of the model, but can induce an extra symmetry in the driven case, as we discuss below.Finally note that driving the system does not break its time-reversal symmetry.We therefore expect a GOE (respectively COE) level statistics in the timeindependent (respectively driven) case [44,73,74,82].Since the COE and GOE ensembles are asymptotically described by the same statistics, we will compare simulations in the driven case to the corresponding GOE statistics.
We choose the parameter set g = Γ × 0.9045, h 2k+1 = 0.809 + 0.9045 × √ 1 − Γ 2 k , Γ = 0.9, τ = π/4, where the k are uniformly distributed random number of zero mean and unit variance.This choice of parameters has been shown [83][84][85] to give good agreement with COE level statistics on the accessible system sizes, for the related model where the longitudinal field is also non-zero on the even sites: h 2k = 0.This is indeed the case for the time-independent system, as can be seen in Fig. 7.
However, the top panel of Fig. 7 shows that there is a dip in the average gap ratio r around J = 1 for the driven system.The numerical estimate for r at this special point appears to coincide with the surmise value for two GOE blocks of equal size given in Table I.The level statistics (bottom panel of Fig. 7) is also compatible with the gap ratio statistics P GOE m=2 (r).Driving the otherwise fully GOE system therefore appears to give rise to a new Z 2 dynamical symmetry at the J = 1 point.
P static (r) P GOE 2 (r) FIG. 7. Top: Average gap ratio as a function of J in the time-independent and driven spin-1/2 chain model Eq.(65).Bottom: Gap ratio distributions at J = 1 for the timeindependent and driven case (points), and comparison with the surmise distribution for m = 1 and m = 2 GOE blocks (full lines).In the time-independent case, average is performed over 2000 disorder realizations (except at the J = 1 point where 5000 realizations are used), and 100 eigenstates around infinite temperature energy E = (Emin + Emax)/2, for the system of size L = 14, while in the driven case, average is performed over 2000 disorder realizations and all eigenstates of the system of size L = 12.
We now rationalize the emergence of this symmetry.We find that the associated conserved operator is Indeed, while this operator acts in general non-trivially on the H z Hamiltonian, we have at the special point J = 1, and thus X commutes with the Floquet operator, up to a global phase factor which can be absorbed in the definition of U F .This indicates the existence of 2 COE blocks of identical sizes in the Floquet Hamiltonian, and hence explains the agreement with the analytically-obtained gap probability distribution P GOE (r).Note that setting h j = 0 on odd (or even) sites is necessary for this extra symmetry to exist.Remark that when open boundary conditions are used, the above commutator Eq. ( 68) becomes proportional to σ z 1 σ z L , a non-trivial boundary term.We can interpret this boundary term as creating a pair of excitations [86].If the model were brought to the MBL regime (e.g. by increasing the strength of the disorder term h j ), these excitations would become localized at both ends of the chain, signalling the topological nature of the observed Z 2 symmetry.However, in that case the level statistics would become Poisson, and we would not be able to detect the symmetry using our RMT approach.Finally, we note that the argument carries over when we add a third contribution H y = L/2−1 k=0 h y σ y 2k+1 .It breaks the time-reversal symmetry, turning the 2-block COE structure of the J = 1 point into a 2-block CUE structure.

D. Anyonic Chain
Our final application deals with chains of interacting anyons, which are exotic particles interpolating between bosonic and fermionic statistics.They are predicted to occur in some two-dimensional systems such as fractional quantum Hall states [87,88], and offer exciting perspectives for topological quantum computation [89].More precisely, we will consider a disordered version of the "golden chain" model of Fibonacci anyons [90].As more technical background is needed to introduce the model and its various representations, we first give a summary of our results.When periodic boundary conditions are imposed on the chain of anyons, there is a non-trivial topological symmetry that decomposes the Hamiltonian into two blocks of unequal size.Resolving this symmetry is not easy, but it can in principle be done at the price of turning the Hamiltonian into a dense matrix, rendering numerical simulations on large systems difficult.Our results instead allow to use a representation simpler for numerics (with real symmetric matrices) which can nevertheless be confronted to RMT predictions, and hence probe ergodic physics.This can be seen in Fig. 8 where the results for r and P (r) allow to characterize the spectral statistics of the model with the two interlaced sectors.At an extra numerical cost, and with the further requirement to study different representations of the model, we can identify the states in each of these two sectors and check that they follow regular single-block GOE statistics (squares and triangles in the top panel of Fig. 8).In the following, we present in detail the different representations of the model of disordered Fibonacci anyons, which allow us to draw these conclusions.
The statistics of anyons [91] are generally encoded in a set of fusion rules analogous to the composition rules for angular momenta, as well as transformation rules relating the different possible ways to fuse together three or more particles (the so-called "F-symbols") [92].In the case of Fibonacci anyons there are only two types of particles, the trivial particle, labeled by 1, and the Fibonacci anyon, labeled by τ .They are characterized by the fusion rule τ × τ = 1 + τ , that is, bringing together two Fibonacci anyons yields either the trivial particle or another Fibonacci anyon.This is analogous to the situation where two spin- 1  2 particles brought together can be decomposed into a spin-0 and a spin-1 particle.In addition the trivial fusion rules 1 × 1 = 1, and 1 × τ = τ hold [90].
Suppose now we have a chain of L indistinguishable Fibonacci anyons.A pair of adjacent anyons may be fused together, yielding either 1 or τ .Performing recursively all possible fusions, we end up with a single anyon, again either 1 or τ .The different ways by which the L particles pair up and fuse to yield a single particle has the structure of a Hilbert space and is called the fusion space.In contrast with the case of spins, this Hilbert space does not have a tensor product structure.In order to construct a basis for this Hilbert space it is convenient to consider the different "fusion paths" which describe the outcome of each fusion, starting from the leftmost pair (particle 1 with particle 2, then the resulting particle with particle 3, and so on).Each fusion path can be written as a sequence |x 1 x 2 . . .x L , where for each i, x i ∈ {1, τ }, and x i+1 is the outcome of the fusion of x i with τ (x 1 being the outcome of the fusion of the first two particles).Since the fusion of 1 with τ always yields τ , no two consecutive 1s are allowed in the sequence of x i .In fact, the basis is given by all strings which do not contain any pair of consecutive 1s.
In the case of periodic boundary conditions, x L+1 = x 1 , and the number of basis states |H| is related to the Fibonacci sequence, |H(L)| = F L−1 + F L+1 , where F L is the Lth Fibonacci number with F 0 = 0 and F 1 = 1.It is well known that the ratio of consecutive Fibonacci numbers goes to lim i→∞ F i+1 /F i = φ with φ = 1+ √ 5 2 the golden mean, hence the name golden chain.A pedagogical introduction to the Hilbert space and Hamiltonian construction of the golden chain model can be found in [93].
Following the seminal work [90], a Hamiltonian can then be constructed by assigning a different energy for each possible outcome of the fusion of two nearest neighbors at sites j and j + 1. Assigning a zero energy to an outcome 1 and −J j to an outcome τ , the Hamiltonian takes the form where Π j,j+1 is the projector into the trivial particle of two τ particles located at sites j and j + 1.This is the analog of the Heisenberg coupling for SU (2) spins 1/2, which assigns a different energy to the fusion channels of pairs of neighboring spins.The projector Π j,j+1 acts on a basis state |x 1 x 2 . . .x L by changing x j to a superposition of 1 and τ in a way depending on x j−1 and x j+1 ; an explicit expression can be found in [90].The coupling constants J j are taken from a random distribution P (J) = −1 J −1+1/ θ(J)θ(1 − J) with θ the Heaviside step function ( ∈ [0, ∞] characterizes the disorder strength and J ∈ [0, 1]).Once again, the main interest for using a disordered coupling constant is to break lattice symmetries in the chain.We use periodic boundary conditions in the following.A practical representation for numerical simulations is to recast the above fusion paths in terms of sequences of heights |h 1 h 2 . . .h L , where h i ∈ {1, 2, 3, 4} and |h i − h i+1 | = 1, through the mapping 1, τ → 1, 3 for i odd, 1, τ → 4, 2 for i even.This defines a "restricted solidon-solid" (RSOS) model, namely the p = 4 case of the A p (also known as SU(2) p−1 ) family, where for generic p the heights are allowed to run between 1 and p [94,95].In this formulation the projectors Π j,j+1 of Eq. ( 69) can be re-expressed in terms of operators e j , whose action is defined as The operators e j form a representation of the Temperley-Lieb (TL) algebra [96], namely e 2 j = √ Qe j , e j e j±1 e j = e j , and e i e j = e j e i for |i − j| ≥ 2, where we have defined √ Q = 2 cos π p+1 .In the present case p = 4, and one indeed checks that Π j,j+1 = 1 √ Q e j .Up to the irrelevant 1/ √ Q proportionality factor, the Hamiltonian Eq. ( 69) is therefore re-expressed in the RSOS representation as A subtlety to keep in mind is that the RSOS formulation acts separately on two equivalent sectors, which correspond to putting even or odd heights on even sites respectively.For periodic boundary conditions, h L+1 = h 1 , each of these sectors has size F L+1 + F L−1 , and yields a copy of the original anyonic chain.The spectrum of the original Hamiltonian Eq. ( 69) is therefore obtained by restricting to a single of these sectors (which is what we do in the following).
The Hamiltonian Eq. ( 71) is real, a reason for which this representation is often used in numerics.We present in Fig. 8 (blue circles) the results for the average gap ratio r and its distribution P (r) for Eq. ( 71), for different chain sizes (and thus Hilbert space sizes |H(L)|) and weak disorder = 0.2, for states located in the middle of the spectrum and corresponding to the sector with even heights on even sites.For this small value of disorder, we do expect a random matrix theory behavior, but the value of r is clearly different from the GOE statistics for a single block.The size of the Hilbert space, which is the sum of two Fibonacci numbers, may suggest the existence of two blocks of different sizes (denoted N 1 and N 2 in the following).A first simple test is to compare the expectation value r RSOS 0.452 to the values obtained for two GOE blocks of different sizes (Fig. 3 in Sec.III D).This leads to a possible value around α = N1 N1+N2 ∈ [0.27, 0.3], corresponding to a size ratio 382.This strongly suggests that the spectrum of the periodic RSOS chain is composed of two independent GOE blocks of size N 1 = F L−1 and N 2 = F L+1 .In the top panel of Fig. 8, we also represent the predicted value r m=2,α=1/(1+φ 2 ) = 0.453186, to which the numerical data indeed appear to tend.This is further confirmed by the numerical distribution of P (r) (bottom panel of Fig. 8) which is an excellent match with the one obtained from the surmise (see Sec. III) of two GOE blocks with ratio φ −2 .
We can in fact trace back this decomposition to the existence of a "hidden" symmetry of a topological nature [90,[97][98][99], namely an operator Y corresponding to an extra τ particle circling around the system, and whose matrix elements in the basis of fusion paths may be written as where the F-symbols F can be found for instance in [90].The operator Y has two distinct eigenvalues 1 2 (1 ± √ 5), and commutes with the Hamiltonian Eq. ( 69), therefore defining two symmetry subspaces.A subtlety arises in the RSOS representation, which as discussed above has two independent sectors and where the action of Y maps one onto the other.We can overcome this difficulty by computing the action of Y 2 , which acts separately in the odd and even sectors: this allows to define in each sector two orthogonal subspaces of dimensions F L+1 and F L−1 respectively, which precisely reproduces the numerical observations made above [100].Now, it is important to remark that the action of Y is highly non-local, and its matrix expression in the RSOS representation is not sparse.Therefore while we know in principle how to decompose the Hamiltonian into two GOE blocks, it is not possible to our best knowledge to do so while keeping it sparse and real.One may ask whether other representations of our model might help with this sectors.The dashed lines represent the results for the (singleblock) GOE distribution of [18] and from the surmise computations in Sec.III for m = 2 GOE blocks with size ratio φ −2 .Data in the Y = 1+ were obtained by considering the rest of the states in the RSOS representation.We focus on midspectrum eigenstates of the RSOS spectrum, obtaining ∼ 300 eigenstates for every disorder realization.Results are averaged over between 300 and 1000 realizations of disorder of strength = 0.2.Bottom: Probability distribution of the gap ratio P (r), as obtained from simulations of a RSOS chain of size L = 22.The solid line is the surmise P GOE m=2,α=1/(1+φ 2 ) (r) obtained from the analytical computations in Sec.III.
problem.There are indeed other ways to represent the TL algebra, from which the spectrum of Eq. ( 69) can be recovered.Below we consider two such representations, the loop representation and the spin chain representation.These representations allow us to tell apart which subspace each eigenvalue corresponds to.
Loop representation In the loop representation [96], the Hilbert space is spanned by the configurations of non-crossing valence bonds between L vertical strands, and the TL generator e i acts by contracting together the strands at site i and i + 1.The composition rules of the TL algebra express the fact that lines can be continuously deformed without crossing, and that closed loops contribute a weight √ Q.From there, one can recover the eigenvalues of the anyon chain corresponding to each symmetry sector by assigning a special weight to noncontractible loops which close around the cylinder [101], respectively 2 cos π 5 and 2 cos 2π 5 , which is nothing but the corresponding eigenvalue of Y .However, the loop model contains significantly more states than the anyonic chain, as the loops carry additional non-local information which is absent in the RSOS representation.This brings several difficulties, the first being that the maximum size L that can be reached using exact diagonalization techniques is lower, the second being that it is not obvious at all how to extract from the loop model spectrum the set of eigenvalues which are present in the RSOS one [102].Furthermore, the loop representation leads to a non-Hermitian matrix representation of the Hamiltonian, which also decreases the efficiency of simulations.
Spin chain representation Another representation is in terms of a spin 1/2 chain, with Hilbert space (C 2 ) ⊗L , on which the TL generators act as [103] Here the matrices σ x,y,z i act as Pauli matrices on the ith spin, and as identity elsewhere, and γ is defined by √ Q = 2 cos γ.The role of the twist parameter ϕ is analog to that of the weights of non-contractible loops in the geometrical representation.More precisely, the Hamiltonian built out of Eq. ( 73) commutes with the global magnetization S z = i σ z i , and the eigenvalues of the RSOS model are recovered in the S z = 0 sector upon setting ϕ = π 5 and ϕ = 2π 5 , respectively.As for the loop case, the S z = 0 sector however contains more states than the RSOS ones, leading to the difficulties mentioned above (see [104][105][106] for other occurrences in related models).Moreover this Hamiltonian is complex in the σ z basis, which also leads to a decreased numerical efficiency.
We use simulations both in the loop and spin chain representations and checked on small systems (up to L = 18) that all states in the RSOS representation can indeed be found in the loop representation (using non-contractible loop weight taking either 2 cos π 5 or 2 cos 2π 5 values) or the spin chain representation (using a twist taking either π/5 or 2π/5 values).The simulations in the loop model with non-contractible loop weight 2 cos π 5 are simpler (as all loops have the same weight) and we could reach larger systems.This allowed us to identify all states in the Y = 1+ sector for chains of size up to L = 24 (see Fig. 8).
Besides allowing to identify this two-block structure in H RSOS chains with periodic boundary conditions, the actual value of r and distribution P (r) for N 1 /N 2 = φ −2 will be useful as a marker of an ETH/ergodic phase when increasing the value of disorder.Indeed, it has been proposed [107] that disordered SU(2) 3 chains could lead to a new form of non-ergodic, critical, phase which behavior is different from a many-body localized phase.This putative new critical phase could be identified by the departure of spectral statistics from the references values displayed above.

V. SUMMARY, RELATION TO PREVIOUS WORKS AND PERSPECTIVES
In summary, we analyzed and computed the statistics of the gap ratio r, an essential tool in diagnosing manybody quantum chaos, when the existence of symmetries results in a block structure of the matrix under consideration.The analytical results we obtain, based on an extension of a seminal work of Rosenzweig and Porter [38], are virtually indistinguishable from numerical simulations on large random matrices.While a closed form can only be obtained in limited cases, our formulation, based on Eqs. ( 27),( 28) and ( 8)), is compact and generic enough to be implemented easily for all cases of interest.Through several examples of applications, we showed the validity and usefulness of our results to identify or probe symmetries in many-body quantum physics.In this final part of this manuscript, we relate our findings to previous works (including a re-interpretation of results available in the literature) and provide leads for possible extensions.

A. Relation to, and re-interpretation of previous results
We now relate our findings to others obtained in studies of spectral statistics in various contexts.Some attempts have been made to count the number of symmetries in chaotic systems [34,35,108].In Appendix A, we provide a comparison between our results and the techniques proposed to detect symmetries in Ref. [108].Our results allow to indirectly discover symmetries in a many-body chaotic system, or to bypass them when they are too complex/costly to implement numerically.There have been several cases of unusual values of P (r) or r reported in previous literature which our work directly elucidates.For instance, it applies to the spectral statistics of the Hamiltonian of the fractional quantum Hall effect when orbital inversion is not resolved in the numerics [109].Our analysis also explains the results obtained on the 2d square lattice quantum Ising model [110] in momentum sectors k = (0, 0) and k = (π, π) where not all symmetries were resolved.The value P (r = 0) 1.4 strongly suggests an unresolved Z 2 symmetry there.Our analysis also accounts for the results in the onedimensional t−t −V clean fermionic model of [111] when the inversion symmetry-breaking field is small, for spectral statistics of the Bose-Hubbard chain [112] when the reflection around the center of the chain is not resolved, as well as of quasiperiodic tilings [113] when phase and parity symmetries are not considered.Another context where our work is relevant is the bosonic SYK model with two-body interactions [26] where the gap ratio distribution (see Fig. 6 in [26]) appears to be close to P GUE m=2 (r), suggesting a two-block GUE structure (for instance due to a particle-hole symmetry), instead of an integrability signature as originally suggested in [26].For some values of the number of Majorana fermions, the bipartite SYK model introduced in Ref. [114] displays the average gap ratio value r 2 that we derive for the GOE ensemble.
In another direction, our analysis could be useful to discover fracton models [115][116][117][118][119] where the Hamiltonian decomposes in several different Krylov independent blocks (and this not necessarily based on an unresolved symmetry), which necessarily implies a non-adherence to the single-block gap ratio statistics [18].A related case is the excellent description of level statistics in an effective quantum ice model [120] with the use of P GOE m=4 (r), accounting for the existing four topological sectors.

B. Perspectives
Our work can be extended in several directions.In a straightforward way, it is possible to extend the analysis to several blocks with different spectral statistics, for instance, the coexistence of GOE and GUE blocks in the same spectrum.This applies to the quantum Hall work of [109] where different momentum sectors have different spectral statistics.Also, it is possible to see the effect of combining integrable and chaotic blocks, in the spirit of the work of [39] on mixed phase spaces.This would apply to models with integrable "Krylov" subspaces coexisting with ergodic blocks [116][117][118], or to the effective model of the MBL transition proposed in [121] with one ergodic block and random independent energies.
Our approach is general enough that it should apply mutatis mutandis to other RMT ensembles or other joint distributions p(s, t).Here we considered the Wigner ensembles with quadratic potential in Eq. ( 1).More generally, β-ensembles with different potentials can be treated in the same way.For instance, β-Laguerre ensembles, with logarithmic potential, are connected with Wishart matrices and are relevant to characterize entanglement spectra (for a review see e.g.[122]).Entanglement spectra can also exhibit block structure inherited from the symmetry of the underlying quantum state from which they are formed.Other natural extensions include replacing Eq. (3), which is our starting point, by the equivalent expression for matrices of larger sizes: indeed [19] obtains, from the exact expression of the joint spacing distribution for 4 × 4 matrices, an expression for P (r) which is more accurate than Eq. ( 7) by an order of magnitude.Another possible direction is to study the non-Hermitian situation [37] with symmetries.
A natural generalization would be to consider higher-order spacing ratios: as was shown numerically in [34,35] higher-order ratios of random matrices with m-fold symmetry can be related with ratios of random matrices with Dyson exponent m, allowing to detect underlying symmetries.An extension to our work could provide analytical grounds to these observations.Finally, it would be interesting to analyze the case of weak symmetry breakings (with small matrix elements between different blocks), using a perturbative approach to estimate P (r) in the same vein as the computation performed for the level spacing distribution in [123].

FIG. 2 .
FIG.2.Distribution of the ratio of consecutive level spacings P (r) for (from top to bottom) GOE, GUE and GSE ensembles and m = 1 to m = 4 blocks (from bottom to top at r = 0).Full lines are the surmises obtained from Sec. III, except the m = 1 case, for which the corresponding surmise Eq. (7) is taken from[18].Points are numerical results from random matrices of size at least Ntot = 2000, averaged over 3.6 × 10 5 histograms.

FIG. 3 .
FIG.3.Top: P GOE 2 (r), for various density ratios α (see text).α φ = 1/(1 + φ 2 ) (in golden color in both plots) is the value that corresponds to the anyonic chain application discussed in Sec.IV D. Bottom: r GOE 2 as a function of α.In both plots, full lines are predictions from the surmises of Sec.III, and points are numerical results obtained on random matrices of size at least Ntot = 2000, averaged over 3.6 × 10 5 realizations.

FIG. 4 . 4 m
FIG. 4. Probability distribution of I 1/4 m for N/m GOE blocks with m = 2 (black) and m = 3 (red).Histograms are obtained from 20000 values, each of which is calculated from 40 realizations of matrices of size N = 180 (solid lines) and from 120 realizations for size N = 48 (dashed lines).Vertical black and red lines indicate the theoretical predictions for m = 2 and 3 respectively; blue line is the mid-value 1 2 (I 1/4 2 + I 1/43 ).

FIG. 6 .
FIG.6.Probability distribution of the gap ratio from the Heisenberg model of Eq. (63) with L = 18 spins, with periodic boundary conditions.Data is averaged over all realizations of the binary transverse field hj = ±1/2 (2 18 in total, 3914 of which are nonequivalent up to symmetries), and over 150 eigenstates around infinite temperature energy E = (Emin + Emax)/2.The red solid line represents the analytical prediction Eq. (64): a linear combination of surmises P GOE

3 PFIG. 8 .
FIG. 8. RSOS model -Top: Average gap ratio r for the RSOS model Eq.(71), as a function of Hilbert space size.Blue circles are the results for the full spectrum of the Hamiltonian (69); squares and triangles are those in the Y = 1± √ 5 2

√ 5 2
sector were obtained by comparing energy spectrum in the loop representation (with noncontractible loop weight 2 cos π/5) and RSOS representation.Data in the Y = 1− √ 5 2