Probing symmetries of quantum many-body systems through gap ratio statistics
PProbing symmetries of quantum many-body systems through gap ratio statistics
Olivier Giraud, ∗ Nicolas Mac´e, † ´Eric Vernier, ‡ and Fabien Alet § Universit Paris-Saclay, CNRS, LPTMS, 91405, Orsay, France Laboratoire de Physique Th´eorique, IRSAMC, Universit´e de Toulouse, CNRS, UPS, France CNRS & LPSM, Universit´e Paris Diderot, place Aur´elie Nemours, 75013 Paris, France (Dated: August 26, 2020)The statistics of gap ratios between consecutive energy levels is a widely used tool, in particular inthe context of many-body physics, to distinguish between chaotic and integrable systems, describedrespectively by Gaussian ensembles of random matrices and Poisson statistics. In this work weextend the study of the gap ratio distribution P ( r ) to the case where discrete symmetries arepresent. This is important, since in certain situations it may be very impractical, or impossible, tosplit the model into symmetry sectors, let alone in cases where the symmetry is not known in thefirst place.Starting from the known expressions for surmises in the Gaussian ensembles, we derive analyticalsurmises for random matrices comprised of several independent blocks. We check our formulaeagainst simulations from large random matrices, showing excellent agreement. We then present alarge set of applications in many-body physics, ranging from quantum clock models and anyonicchains 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 furnishesan efficient way to characterize the number and size of independent symmetry subspaces. We finallydiscuss the relevance of our analysis for existing results in the literature, as well as point out possiblefuture applications and extensions. I. INTRODUCTION
Symmetry considerations are an essential part of aphysicist’s toolbox, with countless applications in allfields of physics, ranging from Noether’s theorem, gaugetheories or the description of phase transitions . An-other frequent tool is the use of simplified models, whichsuccessfully describe the important features of a physi-cal phenomenon without having to deal with all micro-scopic details. In this respect, Random Matrix Theory(RMT), which was first initiated to understand the statis-tical properties of energy levels in complex nuclei , is anextremely successful approach which has also impactedvarious branches in physics . It then comes as no sur-prise that symmetry properties are an integral part ofRMT: one of the best-known examples is the construc-tion of classical Gaussian ensembles from time-reversalsymmetry considerations. Depending on the underlyingsymmetry of the system considered, it is best describedby random matrices belonging to one of the three follow-ing ensembles: Gaussian Orthogonal Ensemble (GOE),Gaussian Unitary Ensemble (GUE) and Gaussian Sym-plectic Ensemble (GSE), whose entries are respectivelyreal, complex or quaternionic random variables. Conve-nient to the description of Floquet operators are the cir-cular ensembles introduced by Dyson : circular Orthog-onal Ensemble (COE), circular Unitary Ensemble (CUE)and circular Symplectic Ensemble (CSE). These ensem-bles have the same asymptotic level spacing distributionsas the Gaussian ensembles .The celebrated conjectures of Berry and Tabor andBohigas, Giannoni and Schmit state that RMT de-scribes the spectral statistics of quantum systems with achaotic semiclassical limit, whereas Poisson statistics pro- vides a description of systems with a classical integrablelimit. These two paradigms serve as reference points tostudy the transitions between localization and ergodic-ity, for instance the Anderson transition as a function ofdisorder . Quite crucially, quantum many-body systems,for which there is in general no semiclassical limit, alsodisplay the same dichotomy: RMT statistics for chaoticsystems and Poisson statistics for quantum integrablesystems, including those showing an emergent integrabil-ity such as many-body localized systems . Numerousexamples illustrate the usefulness of a RMT analysis ofquantum many-body spectra .A universal tool in this respect is the study of the dis-tribution p ( s ) of level spacings, or gaps, defined as the dif-ferences between consecutive energy levels, s i = λ i − λ i − ,assuming that the mean level density is fixed to unity,i.e. (cid:104) s (cid:105) = 1. RMT offers simple, powerful predictions forthe distribution p ( s ) in terms of three different Wignersurmises corresponding to the three Gaussian ensemblesmentioned above . These surmises are obtained by a sim-ple calculation on random 2 × . Normalizing thelevel spacing distribution requires the knowledge of thedensity of states, which is often not analytically avail-able. Numerically, one needs to perform an unfoldingof the spectrum, for which there exists different proce-dures . Unfolding can lead to spurious results , inparticular because of finite-size effects; one may even findinstances where different unfolding procedures leads todifferent physical interpretations of the same data. Formany-body systems, the density of states is genericallyfar from being uniform, which makes the use of the un-folding procedure rather inaccurate. a r X i v : . [ c ond - m a t . d i s - nn ] A ug A very useful alternative to the study of s has beenproposed by Oganesyan and Huse , in terms of the gapratio for three consecutive levels, r i = min( s i ,s i +1 )max( s i ,s i +1 ) . Thekey point is that considering the ratio of gaps ratherthan the gaps themselves suppresses the need to knowor estimate the density of states, and thus avoids thenumerical unfolding step. The probability distribution P ( r ) of gap ratios is thus well-suited to characterize sta-tistical properties of many-body spectra. The Poissonstatistics distribution P Poisson ( r ) = 2 / (1 + r ) can beeasily derived from a Poisson sequence. For the randommatrix spectra, analytical surmises of P GOE ( r ) , P GUE ( r )and P GSE ( r ) have been obtained by Atas et al. fromthe joint eigenvalue distribution of 3 × based on4 × r , in particular its average (cid:104) r (cid:105) andits distribution P ( r ), has become one of the most studiedmetrics in the field of disordered quantum systems. Forinstance, it is often used to characterize the change ofstatistics across a many-body localization (MBL) tran-sition, between an ergodic phase, for which the RMTpredictions for P ( r ) are expected, and a Many-Body Lo-calized phase, which displays emergent integrability andthus P Poisson ( r ) gap ratio statistics . The agreementbetween the RMT-predicted P ( r ) and the numerical es-timate for a given model now routinely diagnoses quan-tum chaotic models. Any discrepancy in the gap ratio asa function of a model parameter is often interpreted as asign of a different physical behavior (see e.g. ). The dis-tribution of gap ratios is also instrumental in analyzingthe symmetry properties of the SYK model and variantsas a function of the number of Majorana fermions . P ( r ) has also been measured experimentally to probe anergodic to MBL transition / crossover . Applications ofthis metrics were also performed in other fields of study,such as in astrophysics or for statistics of the zeros ofthe Riemann zeta function . The computation of thegap ratio statistics has been extended in several ways,such as ratios of gaps for levels with one or more otherlevels in-between, or non-Hermitian matrices .What happens to spectral statistics in the situationwhere symmetries are present in the original Hamilto-nian? In a seminal work , Rosenzweig and Porter com-puted the level spacing distribution P ( s ) of systems withseveral independent random blocks (each being a randommatrix with spacing distribution p ( s )). This situationtypically occurs when a physical system displays discretesymmetries, in which case the number of blocks remainsfinite in the thermodynamic limit. For continuous sym-metries, 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 largeenough systems. In an extension of the original work ,Berry and Robnik considered mixed phase spaces withboth ergodic and integrable blocks . In general, one would be inclined to resolve the under-lying symmetries by treating each block independentlyand performing a block diagonalization. This is not al-ways possible. First, there are cases where a symme-try not previously known or analyzed is discovered fortu-itously (e.g. by monitoring the gap ratio and seeing thatit does not converge to its expected value). Second, insome situations, the block diagonalization is not known,too complex to implement, or cannot be performed to-tally. The later case occurs for instance in systems withnon-Abelian discrete symmetries, where two symmetryoperations that commute with the Hamiltonian do notcommute with each other (see e.g. ). Third, there arecases where the basis transformation leading to a blockstructure in the Hamiltonian is known, but results in aHamiltonian which is more costly to analyze; this is forinstance the case if a sparse Hamiltonian leads to non-sparse blocks, inducing a strong decrease in the perfor-mances of numerical routines (we will present such anexample in Sec. IV D).In this work, we extend the Rosenzweig-Porter anal-ysis to the computation of the gap ratio statistics P ( r )when several independent blocks are present. We do soby calculating the joint gap distribution P ( x, y ) for a ma-trix with several independent random blocks, each beinga random matrix with joint gap distribution p ( s, t ). Weobtain closed expressions for P ( x, y ) and P ( r ) in terms of p ( s, t ) and its primitives. These expressions are valid foran arbitrary number of blocks and an arbitrary distribu-tion p ( s, t ). In the case of Gaussian random matrices, weuse for p ( s, t ) a surmise given by the exact 3 × P ( r ). We note that a recent work provides estimates for P ( r ) and (cid:104) r (cid:105) based on an surmise obtained from explicitanalytical calculations for small-size matrices; however itdoes not take into account all possible level partitions.Our approach is quite different, as we discuss below.Our analytical estimates are virtually indistinguishablefrom numerical simulations on large random matrices.Our results explain several deviations for the distribu-tion P ( r ) or expectation value (cid:104) r (cid:105) observed in the liter-ature, as discussed in Sec. V A. They can also be use-ful in several situations such as those mentioned above(which we illustrate with various applications taken frommany-body physics in Sec. IV), as well as to estimatethe number of effective ergodic blocks in an incompletelythermalized system.The manuscript is organized as follows. We first intro-duce the problem in Sec. II, setting up the notations andsummarizing the useful literature as well as our own re-sults. Sec. III contains the derivation of the generic formof P ( r ) when several independent blocks are present. Wethen present results for the three Gaussian ensembles.Sec. III D compares these analytically obtained resultsto simulations performed on random matrices, showingan excellent agreement. Sec. IV contains several realisticapplications of these results in many-body physics, witha panel of different types of possible symmetries: clocksymmetries, symmetries in disorder realizations, dynam-ical symmetries in Floquet systems, disordered anyonicchains with topological symmetries. In Sec. V, we finallyconclude by first discussing existing examples where ourwork directly applies, and then suggesting some furtherperspectives. II. THE SETTINGA. Random matrix ensembles
Let us first consider the case of a single Gaussian ran-dom matrix H of size N , whose distribution is propor-tional to exp( − tr H ). We denote by λ ≤ . . . ≤ λ N the eigenvalues of such a matrix. The density of eigen-values is given by the Wigner semicircle law ρ ( λ ) = πN √ N − λ . Since there are N ρ ( λ ) δλ levels in aninterval δλ , the corresponding mean level spacing in thevicinity of λ = 0 is ∆ = π/ √ N , which gives a local den-sity 1 / ∆ = √ N /π . The joint distribution of eigenvaluesis P ( λ , . . . , λ N ) = N (cid:89) i Let us now consider ensembles of random matrices ofsize N which can be decomposed into m independentblocks of sizes N , N . . . N m , with (cid:80) i N i = N . Theordered eigenvalues λ < λ < · · · < λ N of such a ma-trix can be obtained by diagonalizing each block sepa-rately and ordering the eigenvalues, so that the spectraof the blocks are interlaced. Let N = ( N , N . . . N m ) bea vector of block sizes. The compound spectrum { λ i , ≤ 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 ˜ r i = s i /s i − , or by the gap ratio distribu-tion P N ( r ), with r i = min( s i /s i − , s i − /s i ) ∈ [0 , P N (˜ r ) = r P N ( r ) holds,which entails that P N ( r ) = 2 P N (˜ r ) θ (1 − ˜ r ) . In thatcase, the distributions of ˜ r and r essentially contain thesame information. As we shall see, this is the case forthe distributions considered in this paper, and therefore,as is often done in numerical simulations, we concentrateon the distribution P N ( r ) with r ∈ [0 , m blocks are independent Gaussian random ma-trices given by the Wigner-Dyson ensembles with index β , then the spectrum of block i , { λ ( i ) q , ≤ q ≤ N i } , ischaracterized by its mean level spacing around λ = 0,given by ∆ i = π/ √ N i , or by its local density ρ i = √ N i /π . The resulting spectrum obtained by the super-position of the m spectra has density ρ = (cid:80) i ρ i . In-troducing the normalized densities µ i = ρ i /ρ , we have µ i = (cid:112) N i /N . For the circular ensembles, where densi-ties are uniform over the unit circle, µ i = N i /N .The Rosenzweig-Porter approach, which gives thenearest-neighbour spacing distribution P ( x ) associatedwith the compound spectrum, consists in assuming thatthe compound spectrum is a superposition of indepen-dent and identically distributed spectra with uniformdensity ρ i and with nearest-neighbour spacing distribu-tion given by the surmise Eq. (2). The computation pro-ceeds by identifying that a gap in the compound spec-trum can originate either from a gap in one of the spectraor from a gap between eigenvalues from two distinct spec-tra. Considering all possibilities and the probabilitiesattached to them leads to the spacing distribution P ( x ).This approach is detailed in Sec. III A. These results wereextended in by considering a mixed phase space, whichamounts to adding Poisson blocks to a chaotic Hamilto-nian. The work derives an explicit formula for P ( x )for a Poisson block and ( N − 1) chaotic blocks with samedensity. C. Summary of our results In this paper we extend the Rosenzweig-Porter ap-proach to derive the joint distribution of consecutivenearest-neighbour spacings P ( x, y ) of a compound spec-trum. It is given by the very compact expressionEq. (27)–(28), for which we give a probabilistic interpre-tation. We then obtain P ( r ) from the analog of Eq. (6),namely P ( r ) = 2 (cid:90) ∞ dx P ( x, rx ) . (8)Applying our expressions to the RMT expressionsEqs. (2) and (3), we obtain a closed general expres-sion for P N ( r ). We then apply this general formula tothe case of identical block sizes N i = N/m , for whichwe use the short notation P m ( r ). Some of these cal-culations result in exact closed (albeit complex) forms,others require a numerical integration. Besides the fulldistribution, we will also consider the average gap ratio (cid:104) r (cid:105) m = (cid:82) rP m ( r ) dr and the limiting value for vanish-ing gap ratio P m (0) = lim r → P m ( r ), as they turn outto be of great practical use to identify the existence of asymmetry ( P (0) = 0 in the no-symmetry case m = 1).Our results are summarised in Table I. In the electronicsupplementary material, we provide a Mathematica note-book allowing to reproduce our calculations. III. ANALYTICAL RESULTSA. Nearest-neighbour spacing distribution P ( x ) Since the function p ( s ) in Eq. (2) corresponds to spec-tra with mean level spacing equal to 1, the spacing dis-tributions of each spectrum are given by the functionEq. (2) rescaled by the mean level spacing, i.e. p ( ρ i s ).We introduce the functions f ( s ) = (cid:90) ∞ da p ( s + a ) (9) m GOE GUE GSE (cid:104) r (cid:105) (cid:104) r (cid:105) m . . . . . . . . . . . . ∞ (Poisson) 0.386294 P m (0) 2 1.40805 1.5228 1.634843 1.71587 1.80758 1.883224 1.83279 1.9023 1.951785 1.88972 1.94334 1.976826 1.92175 1.96413 1.987657 1.94157 1.97582 1.99298 1.95469 1.98292 1.995689 1.96383 1.98748 1.9972410 1.97046 1.99055 1.9981711 1.97541 1.99269 1.9987512 1.97922 1.99423 1.99912 . . . . . . . . . . . . ∞ (Poisson) 2TABLE I. Values of averages (cid:104) r (cid:105) and probability at r = 0for m blocks, obtained from the surmise approach in Sec. III.The value for m = 1 is taken from . Values for (cid:104) r (cid:105) m obtainedfrom numerical simulations of random matrices are presentedin Tab. II in Sec. III D. and g ( s ) = (cid:90) ∞ da (cid:90) ∞ db p ( s + a + b ) . (10)The function f ( s ) gives the probability to have λ i +1 ≥ s knowing that λ i = 0. The function g ( s ) gives the proba-bility to have λ i +1 ≥ s knowing that λ i ≤ 0, that is, theprobability to have a spacing at least s . These probabil-ities are related through the identities g (cid:48) = − f, g (cid:48)(cid:48) = p. (11)Introducing the rescaled spacing x = ρs we have fromEq. (11) g ( ρ i s ) = g ( µ i x ) ,µ i f ( ρ i s ) = − ∂ x g ( µ i x ) ,µ i p ( ρ i s ) = ∂ x g ( µ i x ) . (12)Spacings arise as empty intervals ] λ ( i ) q , λ ( j ) q (cid:48) [ of length s with i, j = 1 , ..., m . We have to consider the two pos-sibilities i = j or i (cid:54) = j , and calculate the probabilitydensities associated with each configuration.Let us first consider the case where m = 2 (we willlater generalize this analysis to more blocks). We haveto consider the two following cases:1. Configurations giving rise to an empty interval oftype ] λ ( i ) q , λ ( i ) q +1 [, which are are such that λ ( j ) q (cid:48) <λ ( i ) q < λ ( i ) q +1 < λ ( j ) q (cid:48) +1 for some q (cid:48) . The prob-ability of such a configuration for i is given by p ( ρ i s ), while the probability for j is g ( ρ j s ) since λ ( j ) q (cid:48) and λ ( j ) q (cid:48) +1 can be anywhere outside ] λ ( i ) q , λ ( i ) q +1 [.Taking into account the probability µ i to have alevel i at both ends of the interval, we get forthe configuration ] λ ( i ) q , λ ( i ) q +1 [ a probability density µ i p ( ρ i s ) g ( ρ j s ). Using Eq. (12) we can rewrite it as[ ∂ x g ( µ i x )] g ( µ j x ).2. Configurations giving rise to an empty interval oftype ] λ ( i ) q , λ ( j ) q (cid:48) [, which are such that λ ( j ) q (cid:48) − < λ ( i ) q <λ ( j ) q (cid:48) < λ ( i ) q +1 . The probability of such a spacingfor i is given by f ( ρ i s ) since λ ( i ) q +1 can be any-where in ] λ ( j ) q (cid:48) , ∞ [, while the probability for j is f ( ρ j s ) since λ ( j ) q (cid:48) − 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 thatthe probability density of configuration ] λ ( i ) q , λ ( j ) q (cid:48) [ is µ i µ j f ( ρ i s ) f ( ρ j s ). Using Eq. (12) it can be rewrit-ten as [ ∂ x g ( µ i x )] [ ∂ x g ( µ j x )].Summing these probabilities over i, j = 1 , P ( x ) of the mixed levels P ( x ) = g ( µ x ) ∂ x g ( µ x ) + g ( µ x ) ∂ x g ( µ x )+ 2 [ ∂ x g ( µ x )] [ ∂ x g ( µ x )] (13)or equivalently P ( x ) = ∂ x G ( x ) , G ( x ) = m (cid:89) i =1 g ( µ i x ) . (14)The above reasoning proves Eq. (14) for m = 2. Infact, given the final expression, we can come up with amuch shorter proof of the validity of Eq. (14) for arbi-trary m . Indeed, the probability of finding an intervalof a given length between two consecutive eigenvalues isthe second derivative of the probability of finding an in-terval larger or equal to it with no eigenvalue in it (seeEqs. (10) and (11)). Therefore P ( x ) must be the secondderivative of the probability of finding an empty inter-val larger or equal to x . The probability that no level oftype i occurs in an interval of size x is g ( µ i x ). Since levelsfrom different sequences are independent, the probabilitythat no level of any type occur in an interval of size x issimply the product of all g ( µ i x ), which directly entailsEq. (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 jointdistribution of two consecutive spacings. Our aim is to obtain the joint distribution P ( x, y ) in terms of the dis-tribution p ( s, t ) for a single spectrum.Starting with the function p ( s, t ) in Eq. (3), we intro-duce the functionˆ p ( s ) = (cid:90) ∞ da p ( s, a ) , (15)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 byintegrating over the first variable. Note that this expres-sion differs from the one in Eq. (2), which correspondsto the result for 2 × p ( s, t ) was obtainedfor 3 × f and g can then be de-fined from ˆ p as in Eq. (9)–(10). In terms of p ( s, t ), theirexplicit form is f ( s ) = (cid:90) ∞ da (cid:90) ∞ db p ( s + a, b ) (16)and g ( s ) = (cid:90) ∞ da (cid:90) ∞ db (cid:90) ∞ dc p ( s + a + b, c ) . (17)We also define the two-variable functions e ( s, t ) = (cid:90) ∞ da p ( s + a, t ) , e ( s, t ) = (cid:90) ∞ da p ( s, t + a ) , (18)and h ( s, t ) = (cid:90) ∞ da (cid:90) ∞ db p ( s + a, t + b ) . (19)These functions are related through the identities ∂ s ∂ t h = p, ∂ s h = − e , ∂ t h = − e , g (cid:48) = − f, g (cid:48)(cid:48) = ˆ p. (20)The analogs of Eqs. (12) are g ( ρ i s ) = g ( µ i x ) , µ i f ( ρ i s ) = − ∂ x g ( µ i x ) , (21) µ i p ( ρ i s, ρ i t ) = ∂ x ∂ y h ( µ i x, µ i y ) (22) ∂ x h ( µ i x, µ i y ) = − µ i e ( ρ i s, ρ i t ) (23) ∂ y h ( µ i x, µ i y ) = − µ i e ( ρ i s, ρ i t ) . (24)A sequence of two consecutive spacings arises as a se-quence of two empty intervals ] λ ( i ) q , λ ( j ) q (cid:48) [ and ] λ ( j ) q (cid:48) , λ ( k ) q (cid:48)(cid:48) , [with i, j, k = 1 , ..., m . We have to consider all possi-bilities for i, j, k and calculate the probability densitiesassociated with each configuration.Once again, let us consider the simplest case m = 3,that we will later generalize. We only need to examinefour cases, corresponding to patterns iii , iij , iji and ijk and depicted in Fig. 1 :1. Configurations λ ( i ) q < λ ( i ) q +1 < λ ( i ) q +2 , which arisewhenever λ ( j ) q (cid:48) < λ ( i ) q < λ ( i ) q +1 < λ ( i ) q +2 < λ ( j ) q (cid:48) +1 forsome q (cid:48) and λ ( k ) q (cid:48)(cid:48) < λ ( i ) q < λ ( i ) q +1 < λ ( i ) q +2 < λ ( k ) q (cid:48)(cid:48) +1 for λ ( i ) q λ ( j ) q (cid:48) λ ( k ) q (cid:48)(cid:48) FIG. 1. The four configurations of consecutive spacings con-sidered in Sec. III B. Different colours correspond to distinctspectra. The three central levels are the ones from which theratios are calculated; the outer levels are at the same heightto indicate that their relative position is irrelevant. some q (cid:48)(cid:48) . The probability of such a configurationfor i is given by p ( ρ i s, ρ i t ); the probability for j is g ( ρ j ( s + t )) since λ ( j ) q (cid:48) and λ ( j ) q (cid:48) +1 can be anywhereoutside ] λ ( i ) q , λ ( i ) q +2 [; and the same goes for k . Takinginto account the probability µ i to have three levels i at the ends of the intervals, we get for this config-uration a probability density µ i p ( ρ i s, ρ i t ) g ( ρ j ( s + t )) g ( ρ k ( s + t )). Using Eqs. (21)–(24) we can rewriteit as µ i [ ∂ x ∂ y h ( µ i x, µ i y )] g ( µ j ( x + y )) g ( µ k ( x + y )).2. Configurations λ ( i ) q < λ ( i ) q +1 < λ ( j ) q (cid:48) , which arisewhenever λ ( j ) q (cid:48) − < λ ( i ) q < λ ( i ) q +1 < λ ( j ) q (cid:48) < λ ( i ) q +2 and λ ( k ) q (cid:48)(cid:48) < λ ( i ) q < λ ( i ) q +1 < λ ( j ) q (cid:48) < λ ( k ) q (cid:48)(cid:48) +1 for some q (cid:48)(cid:48) .The probability for i is e ( ρ i s, ρ i t ), the probabilityfor j is f ( ρ j ( s + t )), while the probability for k is g ( ρ k ( s + t )). We get for this configuration a prob-ability density µ i µ j e ( ρ i s, ρ i t ) f ( ρ j ( s + t )) g ( ρ k ( s + t )) = µ i [ ∂ x h ( µ i x, µ i y )][ ∂ y g ( µ j ( x + y )] g ( µ k ( x + y )).We used the fact that ∂ y g ( x + y ) = g (cid:48) ( x + y ).3. Configurations λ ( i ) q < λ ( j ) q (cid:48) < λ ( i ) q +1 which arisewhenever λ ( j ) q (cid:48) − < λ ( i ) q < λ ( j ) q (cid:48) < λ ( i ) q +1 < λ ( j ) q (cid:48) +1 and λ ( k ) q (cid:48)(cid:48) < λ ( i ) q < λ ( j ) q (cid:48) < λ ( i ) q +1 < λ ( k ) q (cid:48)(cid:48) +1 for some q (cid:48)(cid:48) .The probability for i is ˆ p ( ρ i ( s + t )), the probabil-ity for j is h ( ρ j s, ρ j t ), and the probability for k is g ( ρ k ( s + t )). We get for this configuration a prob-ability density µ i µ j ˆ p ( ρ i ( s + t )) h ( ρ j s, ρ j t ) g ( ρ k ( s + t )) = µ j [ ∂ x ∂ y g ( µ i ( x + y ))] h ( µ j x, µ j y ) g ( µ k ( x + y )).We used the fact that ∂ x ∂ y g ( x + y ) = g (cid:48)(cid:48) ( x + y ).4. Finally, configurations λ ( i ) q < λ ( j ) q (cid:48) < λ ( k ) q (cid:48)(cid:48) , whicharise whenever λ ( j ) q (cid:48) − < λ ( i ) q < λ ( j ) q (cid:48) < λ ( k ) q (cid:48)(cid:48) < λ ( j ) q (cid:48) +1 and λ ( k ) q (cid:48)(cid:48) − < λ ( i ) q < λ ( j ) q (cid:48) < λ ( k ) q (cid:48)(cid:48) < λ ( i ) q +1 . The prob-ability 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 den-sity µ i µ j µ k f ( ρ i ( s + t )) h ( ρ j s, ρ j t ) f ( ρ k ( s + t )) = µ j [ ∂ x g ( µ i ( x + y ))] h ( µ j x, µ j y )[ ∂ y g ( µ k ( x + y ))].We can now sum all contributions over i, j, k = 1 , , P ( x, y ) = ∂ x ∂ y [ µ h ( µ x, µ y ) g ( µ ( x + y )) g ( µ ( x + y ))+ µ g ( µ ( x + y )) h ( µ x, µ y ) g ( µ ( x + y ))+ µ g ( µ ( x + y )) g ( µ ( x + y )) h ( µ x, µ y )] . (25)Equation (25) has a simple probabilistic interpretation.Let us define a function H as H ( x, y ) = (cid:90) ∞ da (cid:90) ∞ db P ( x + a, y + b ) , (26)by analogy with Eq. (19). Thus H ( x, y ) gives the prob-ability of having a triplet ( λ q − , λ q , λ q +1 ) of levels ofthe mixed spectrum such that λ q − < λ q − x and λ q + y < λ q +1 . That is, H ( x, y ) is the probability thatsome level λ ( i ) q is such that all other levels λ ( j ) q (cid:48) verifyeither λ ( j ) q (cid:48) < λ ( i ) q − x or λ ( i ) q + x < λ ( j ) q (cid:48) (including thecase i = j , in which case of course q (cid:48) (cid:54) = q ). At fixed i ,the probability of such a configuration is h ( µ i x, µ i y ) forspectrum i , and g ( µ j ( x + y )) for all j (cid:54) = i . Summing overall i (and taking into account the probability µ i to havea level i in the middle), we get for H ( x, y ) the expressionunder the derivation symbols in Eq. (25).In fact, this reasoning provides a proof of the generalcase with arbitrary number m of spectra. We thus havein the general case P ( x, y ) = ∂ x ∂ y H ( x, y ) (27)with H ( x, y ) = m (cid:88) i =1 µ i h ( µ i x, µ i y ) (cid:89) j (cid:54) = i g ( µ j ( x + y )) . (28)One can check, using Eqs. (20)–(24), that P ( x, y ) inEq. (27) is properly normalized to 1. C. Ratio distribution P ( r ) In order to obtain P ( r ), one first needs to evaluatefunctions g and h to obtain H ( x, y ) using Eq. (28), thentake its derivative with respect to x and y , and finallyperform the integration in Eq. (8). In the GUE and GSEcase, there is no closed-form expressions for the function h , so that we are left with a double integral (one in thedefinition of h , one corresponding to the final integrationin Eq. (8)). In the GOE case however, we obtain explicitexpressions for g and h , as we will show below, and thuswe get a closed-form expression for P ( x, y ). The remain-ing integral Eq. (8) is doable analytically only in the caseof a mixture of m = 2 spectra. 1. GOE case In the GOE case, the joint distribution p ( s, t ) reads p ( s, t ) = 3 π st ( s + t ) e − π ( s + st + t ) . (29)Starting from Eq. (29) for p ( s, t ) and calculating explic-itly the functions g and h given in Eqs. (17) and (19) weget g ( s ) = U ( s ) − s U ( s ) − s U ( s ) (30)with U ( s ) = e − π s ,U ( s ) = Erfc (cid:18) s √ π (cid:19) ,U ( s ) = e − π s Erfc (cid:18) s √ π (cid:19) , (31)and h ( s, t ) = 9( s + t )4 π V ( s, t )+ 8 π − s π V ( s, t ) + 8 π − t π V ( s, t ) (32)with V ( s, t ) = e − π ( s + st + t ) ,V ( s, t ) = e − s π Erfc (cid:18) s + 2 t )4 √ π (cid:19) ,V ( s, t ) = e − t π Erfc (cid:18) s + t )4 √ π (cid:19) . (33)One can then rewrite Eq. (28) as H ( x, y ) = m (cid:88) i =1 3 (cid:88) a =1 · · · (cid:88) a m =1 H ( i ) a ...a m V a i ( µ i x, µ i y ) (cid:89) j (cid:54) = i U a j ( µ j ( x + y ))(34)with H ( i ) a ...a m some polynomials of x , y and the µ i , whichcan be obtained explicitly from Eqs. (30) and (32).Functions U i have the property that they transforminto each other under derivation. Namely, the derivativeof any function (cid:80) i c i U i (with c i polynomial in s ) is ofthe form (cid:80) i ˜ c i U i (with ˜ c i polynomial in s ). The sameproperty holds for the V i upon derivation with respectto s or t . Therefore, using Eq. (27) and the expansionEq. (34), we obtain P ( x, y ) = m (cid:88) i =1 3 (cid:88) a =1 · · · (cid:88) a m =1 P ( i ) a ...a m V a i ( µ i x, µ i y ) (cid:89) j (cid:54) = i U a j ( µ j ( x + y ))(35)where P ( i ) a ...a m are polynomials of x , y and the µ i . Giventhe definition of the functions U i and V i , Eq. (8) can beexpanded as a linear combination (with real coefficientsdependent on the µ i ) of integrals of the form (cid:90) ∞ dx x k e − λx m (cid:89) i =1 Erfc( a i x ) , (36) with λ and the a i depending on the µ i and on r (andpossibly a i = 0). It appears that in general such anintegral does not have a closed form. However in thecase m = 2 we have the identity (cid:90) ∞ dx xe − λx Erfc( ux ) Erfc( vx ) =12 λ − u tan − (cid:16) √ λ + u v (cid:17) πλ √ λ + u − v tan − (cid:16) √ λ + v u (cid:17) πλ √ λ + v , (37)from which one can deduce Eq. (36) for all odd values of k by deriving with respect to λ , and for all even valuesof k by first integrating by parts and then deriving withrespect to λ . This yields a (rather lengthy) closed-formexpression for P ( r ) in the case m = 2 (which is given infull in the electronic Supplementary Material). To givean idea of this expression, we only give P (0) in the caseof two blocks of the same size: P (0) = 1168 (cid:16) − √ √ π +14 √ − (cid:18) √ (cid:19) − √ − (cid:18) √ (cid:19)(cid:19) (cid:39) . . (38)For m blocks of the same size, we get P (0) = 1 . m = 3 and P (0) = 1 . m = 4, as reported inTab. I.In practice, the fastest way of obtaining P ( r ) for GOEin 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 nu-merically. The Mathematica notebook in the electronicSupplementary Material implements the two possibilitiesto obtain P ( r ). From the exact equation, numerical in-tegration over [0 , 1] yields the mean ratio. For instancefor m blocks of equal size, we get (cid:104) r (cid:105) GOE ,m blocks = 0 . , . , . m = 2 , 2. GUE case In the GUE case, the joint distribution p ( s, t ) reads p ( s, t ) = 3 √ π s t ( s + t ) e − π ( s + st + t ) . (40)The calculation of g and h for GOE was made possible bythe fact that either s or t is of degree 1 in the polynomialin front of the exponential in Eq. (29). This is no longerthe case for GUE and GSE. However, g as well as the firstderivative of h can be obtained analytically. We find g ( x ) = − √ x π e − x π − Erfc (cid:32) (cid:114) π x (cid:33) + e − x π (cid:32) (cid:0) x + 128 πx (cid:1) π + 2 (cid:33) Erfc (cid:32) (cid:114) π x (cid:33) − x (cid:18) Erfc (cid:18) x √ π (cid:19) − T (cid:18) x √ π , √ (cid:19)(cid:19) (41)where T ( x, a ) is Owen’s T -function, defined as T ( x, a ) = 12 π (cid:90) a dt 11 + t e − x (1+ t ) / . (42)For h we have h ( x, y ) = (cid:82) ∞ e ( x + a, y ) da , with e ( x, y ) = (cid:90) ∞ db p ( x, y + b )= 3 √ x e − ( x xy + y ) π ( x + 2 y )2 π × (cid:0) π − (cid:0) x − xy − y (cid:1)(cid:1) + 3 e − x π x Erfc (cid:16) (cid:113) π ( x + 2 y ) (cid:17) π × (cid:0) x − πx + 16384 π (cid:1) (43)and e ( x, y ) = e ( y, x ). Using the identities in Eq. (20),the expression of P ( x, y ) reads P ( x, y ) = m (cid:88) i =1 µ i h ( µ i x, µ i y ) ∂ x ∂ y (cid:89) j (cid:54) = i g ( µ j ( x + y )) − m (cid:88) i =1 µ i e ( µ i x, µ i y ) ∂ y (cid:89) j (cid:54) = i g ( µ j ( x + y )) − m (cid:88) i =1 µ i e ( µ i x, µ i y ) ∂ x (cid:89) j (cid:54) = i g ( µ j ( x + y )) + m (cid:88) i =1 µ i p ( µ i x, µ i y ) (cid:89) j (cid:54) = i g ( µ j ( x + y )) , (44)in which only the first sum involves the function h forwhich no closed form is available. The i th term in thatsum is µ i h ( µ i x, µ i y ) q i ( x + y ), with q i an explicitly knownfunction involving only products of derivatives of g . Forthis term the integral Eq. (8) yields a contribution2 µ i (cid:90) ∞ dx x h ( µ i x, µ i rx ) q i ((1 + r ) x )= 2 µ i (cid:90) ∞ dx (cid:90) ∞ da x e ( µ i x + a, µ i rx ) q i ((1 + r ) x ) . (45) 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: P (0) = 1 . , . , . (cid:104) r (cid:105) GUE ,m blocks = 0 . , . , . m = 2 , , 3. GSE case For GSE the joint distribution p ( s, t ) reads p ( s, t ) = 3 √ π s t ( s + t ) e − π ( s + st + t ) . (48)The function g reads g ( x ) = − x (cid:18) Erfc (cid:18) x √ π (cid:19) − T (cid:18) x √ π , √ (cid:19)(cid:19) − Erfc (cid:32) (cid:114) π x (cid:33) + e − x π Erfc (cid:16) (cid:113) π x (cid:17) Q ( x ) − √ e − x π x π Q ( x ) (49)with Q , Q the polynomials Q ( x ) = 523347633027360537213511521 x π + 8862938119652501095929 x π + 16677181699666569 x π + 2792914305201 x π + 14703201 x π + 20736 (50)and Q ( x ) = 16677181699666569 x π + 94143178827 x π + 531441 x π − . (51)It can be checked that the second derivative of g is indeedthe marginal probability ˆ p , and that g (0) = 1. Similarlyas in the GUE case, h is not calculable in closed form,but we have h ( x, y ) = (cid:82) ∞ e ( x + a, y ) da with e ( x, y ) = 3 x π e − ( x xy + y ) π (cid:32) R ( x, y )+ R ( x, y ) e x +2 y )2102400 π Erf (cid:32) (cid:114) π ( x + 2 y ) (cid:33)(cid:33) , (52)where R ( x, y ) and R ( x, y ) are polynomials given by R ( x, y ) = 77760 √ x + 2 y ) × (cid:18) − π (cid:0) x − xy − y (cid:1) +535570083993600 π (cid:0) x − x y + 88 x y + 224 xy + 112 y (cid:1) − (cid:0) x − x y + 12 x y − x y − x y − xy − y (cid:1) +4697620480000000 π (cid:19) (53)and R ( x, y ) = 328256967394537077627 x − πx + 493581389408501760000 π x − π x + 240518168576000000000 π . (54)The computation is then the same as for GUE. UsingEqs. (44) and (45) we get, for m = 2 , , P (0) = 1 . , . , . (cid:104) r (cid:105) GSE ,m blocks = 0 . , . , . . (56)We finally mention that the Mathematica notebook inthe electronic Supplementary Material allows to repro-duce these computations, as well as to consider differentcases (other values of m , unequal sizes of the m blocks). 4. Poisson ( m → ∞ ) limit One can easily check that in the case of a mixtureof m spectra of the same size, one recovers the Poissondistribution in the m → ∞ limit. Indeed, in that casethe function H reads H ( x, y ) = h (cid:16) xm , ym (cid:17) g (cid:18) (1 + r ) xm (cid:19) m − (57)and thus P ( r ) = 2 (cid:90) ∞ dx x g m − (cid:20) ∂ x ∂ y hm + m − m ( ∂ x h + ∂ y h ) g (cid:48) g + h m − m (cid:18) ( m − g (cid:48) g + g (cid:48)(cid:48) g (cid:19)(cid:21) , (58)with functions h and g evaluated at (cid:0) xm , rxm (cid:1) and (1+ r ) xm ,respectively. In the limit m → ∞ , these argumentsgo to 0. From the explicit expressions for g and h , the only term in the square brackets that survives is h ( m − m − m g (cid:48) g , which goes to 1. Using the fact that g ( x ) = 1 − x + O ( x ) close to 0, we get P ( r ) (cid:39) (cid:90) ∞ dx x (cid:18) − (1 + r ) xm (cid:19) m − → m →∞ (cid:90) ∞ dx x e − (1+ r ) x = 2(1 + r ) , (59)which is indeed the Poisson result. D. Comparison with numerics P ( r ) GOE m = 1 m = 2 m = 3 m = 4 P ( r ) GUE0 . . . . . . r P ( r ) GSE FIG. 2. Distribution of the ratio of consecutive level spacings P ( r ) for (from top to bottom) GOE, GUE and GSE ensemblesand 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 . Points are numerical results from randommatrices of size at least N tot = 2000, averaged over 3 . × histograms. We now compare the gap ratio distribution P m ( r ) andthe mean value (cid:104) r (cid:105) m obtained through the analytical ap-proach of Sec. III to direct numerical computations on0large random matrices. Numerical RMT spectra are com-puted using the matrix models of , based on tridiag-onal matrices. These models are numerically faster todiagonalize, but otherwise equivalent to the dense ones.By construction, the spectrum of a single RMT block oflinear size N has support in [ − √ N , √ N ]. The sup-ports of a collection of blocks therefore overlap in the[ − √ N min , √ N min ] region, where N min = min j N j is thelinear size of the smallest block in the collection. In orderto avoid boundary effects, we restrict the numerical com-putation of the level statistics to the central quarter ofthis overlap region. The normalized densities µ i in thatregion are given by µ i = (cid:113) N i N tot with N tot = (cid:80) j N j thetotal linear size. In all numerical computations presentedhere, N tot is at least 2 × , and 3 . × samples areused.Figure 2 displays the results of this comparison for m = 2 , , for m = 1)and all Gaussian ensembles. The comparison is excel-lent and within the scale of this figure, there is no visibledifference between the analytically obtained P m ( r ) andthe numerical results P num m ( r ). More precisely, we foundthat the relative error | P m ( r ) − P num m ( r ) | /P m ( r ) is alwaysless than 0 . 01. This translates into also an almost perfectagreement (with no difference within error bars) between (cid:104) r (cid:105) m (from Tab. I) and the numerical estimates reportedin Tab. II. m GOE GUE GSE (cid:104) r (cid:105) (cid:104) r (cid:105) m (cid:104) r (cid:105) for m blocks, as obtainedfrom simulations on random matrices. The value for m = 1 istaken from . Notice the excellent agreement with the surmiseresults reported in Tab. I. We also compare analytical and numerical results forthe specific case of m = 2 GOE blocks of different sizes.In Fig. 3 we present results for the probability distribu-tion P ( r ) for different values of α = µ µ + µ (top panel) aswell as for the expectation value (cid:104) r (cid:105) ,α as a function of α (bottom panel). Here again the agreement between ana-lytical 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 α ≥ . 2. For α < . 2, the relative error increases, but remains below0 . 05. The seemingly crossing point in Fig. 3 is in fact nota crossing point, as one can convince oneself by using theexact expression for P GOE2 ( r ) and calculating its valuewith enough precision in the vicinity of that point.We conclude this comparison section by discussing ,which provides an analytical estimate for P ( r ) for m = 2,derived from a 4 × m = 2 blocks. This does notcontain all possible patterns considered in Sec. III B. The . . . . . . r . . . . P ( r ) α = 0 . α = 0 . α = 0 . α = α φ α = 0 . α = 0 . α = 0 . . . . . . . α . . . . . h r i G O E h r i GOE2 ,α h r i RMT2 ,α FIG. 3. Top: P GOE2 ( r ), for various density ratios α (see text). α φ = 1 / (1 + φ ) (in golden color in both plots) is the valuethat corresponds to the anyonic chain application discussed inSec. IV D. Bottom: (cid:104) r (cid:105) GOE2 as a function of α . In both plots,full lines are predictions from the surmises of Sec. III, andpoints are numerical results obtained on random matrices ofsize at least N tot = 2000, averaged over 3 . × realizations. estimated value for (cid:104) r (cid:105) obtained from this approach isapproximately close to the one presented in Table I forthe GOE, GUE but strongly differs for the GSE, whileour results in this latter case agree with the numericalestimates in Table II. IV. ILLUSTRATIONS IN QUANTUMMANY-BODY PHYSICS We now illustrate the usefulness of the above resultsby comparing them with simulations on realistic spec-tra obtained from quantum many-body problems. Mostof our examples are taken from one-dimensional latticemodels, mostly for computational convenience. In thefollowing, the lattice will thus be a one-dimensional chainwith L sites. Except otherwise mentioned, we will explic-itly break translation symmetry, as well as possibly otherlattice symmetries (such as reflection around the centerof the chain) to concentrate on the existence of a fewblocks. The existence of translation symmetry would re-sult in the existence of L blocks (labeled by the L recipro-1cal wave-numbers), which would result, as discussed ear-lier, in an (effective) Poisson distribution for level spac-ings and gap ratios in the thermodynamic limit. Thetranslation symmetry will be broken by using disordercharacterized by a disorder strength (cid:15) . In all the simula-tions presented below, we take (cid:15) not too small (in orderto avoid the proximity to the translation-invariant case,which would cause stronger finite-size effects) as well asnot too large, to avoid for instance a possible many-bodylocalized phase (which would also result in Poisson spec-tral statistics). In all Hamiltonian systems we examine,we consider mid-spectrum eigenstates, obtained either byfull diagonalization (for the smaller Hilbert space sizes)or by the shift-invert subset diagonalization method forlarger systems. A. Quantum clock models The first example deals with Q -states quantum clockmodels, which are natural Z Q -symmetric generaliza-tions of the Ising quantum chain with Q -states quantum“spins” on each site . These exhibit a rich groundstate phase diagram including ordered and disorderedphases as well as critical lines, and have attracted alot of attention in the recent years due to their rela-tion with parafermions, a Z Q generalization of Majoranafermions , as well as with topological phases . Theyare furthermore related to cornerstone models of statisti-cal mechanics, including the Potts model (where the Z Q symmetry is promoted to a larger, S Q symmetry) andthe chiral Potts model .On each site, we define a spin taking Q possibles values(0 . . . Q − σ and τ , whichgeneralize the Pauli matrices σ z and σ x of the Ising chain: σ measures the orientation of the spin, while τ rotates itby one unit “around the clock”, and as a result thesefulfill the following algebraic rules: σ Q = τ Q = 1, σ † = σ Q − , τ † = τ Q − and στ = ωτ σ with ω = exp(2 iπ/Q ),a Q th root of unity.Simple matrix representations are obtained in the basiswhere σ is diagonal (the “ { σ } -basis”): σ = ω . . . ω Q − , τ = . (60)In the basis where τ instead is diagonal (the “ { τ } -basis”),the matrices are exchanged.The standard Hamiltonian for quantum clock mod-els is written as a linear combination of ( τ j ) a , a =1 , . . . , Q − j and exchange terms ( σ † j σ j +1 ) a , a = 1 , . . . , Q − 1. It is invariant under a Z Q “clock” sym-metry σ j → ωσ j , and the associated conserved charge Z = (cid:81) j τ j has eigenvalues 1 , ω, . . . ω Q − . For Q ≥ τ j → τ † j , σ j → σ † j , andtime-reversal, which is anti-unitary (and therefore sendsany constant to its complex conjugate) and sends σ j to σ † j while leaving τ j invariant.The model we consider in the following breaks all sym-metries, but Z Q : H Q = − (cid:88) j J j σ † j σ j +1 + Γ τ j + ig ( τ j − τ † j ) σ j σ † j +1 + h.c., (61)where the sum runs over the L sites of the 1d lat-tice. For practical computations we restrict ourselves to Q = 2 , , 4. The coupling constants J j are independentrandom numbers uniformly taken from a box distribution[ J − (cid:15), J + (cid:15) ]. Since they are a priori different on eachsite, they break invariance under translation or spatialreflection. The last term breaks both time-reversal andcharge conjugation symmetry (this breaking could alsohave been achieved by perturbing with the U (1) charge S z introduced in ).We first consider results of simulations performed inthe { σ } basis, with a full Hilbert space of size Q L . Inthe top panel of Fig. 4 we present the average gap ra-tio for different chain sizes. We clearly observe gap ra-tios which do not tend to their GOE value (cid:104) r (cid:105) GOE , butrather to their (cid:104) r (cid:105) GOE m = Q value as the size of the Hilbertspace is increased. This is expected, as the Hamiltonianpossesses Q sectors of identical size Q L − labeled withthe different eigenvalues of the charge Z . Note howeverthat working in the { σ } basis does not allow to simplyconstruct the Hamiltonian blocks, as Z is off-diagonal inthat case. Furthermore, the Hamiltonian is complex inthis basis, and without any further indication on the ex-istence of the Z Q symmetry, it would not be clear whyGUE statistics should not show up.When switching to the { τ } basis, Z is now diagonaland the blocks are easily constructed. Furthermore theHamiltonian becomes real. Computing the average gapratio in each block leads to an asymptotic (cid:104) r (cid:105) GOE valuefor each block, showing that each block is indeed inde-pendent and no further symmetry has been missed.We further confirm these results by showing the fulldistribution P ( r ) in the bottom panel of Fig. 4 for Q = 2 , , 4. When the full spectrum is taken, the distri-bution obtained numerically for the largest system sizeis in excellent agreement with the surmises P GOE m = Q ( r ) ob-tained from Sec. III.Note the importance of the time-reversal symmetrybreaking term g (cid:54) = 0 in Eq. (61) in this analysis. Inthe presence of time-reversal symmetry (at g = 0), theblocks with Z and Z ∗ are identical, leading to exact de-generacies. These additional values at r = 0 would resultin effective values of (cid:104) r (cid:105) lower than their Poisson valuesfor finite-size systems.2 |H| . . . . . . . . . h r i h r i GOE h r i GOE2 h r i GOE3 h r i GOE4 h r i Poisson Q = 2 Q = 3 Q = 4 . . . . . . r . . . . . . . . . P ( r ) P ( r ) GOE m =2 P ( r ) GOE m =3 P ( r ) GOE m =4 Q = 2 Q = 3 Q = 4 FIG. 4. Q-states clock model — Top: Average gap ratio (cid:104) r (cid:105) for the model of Eq. (61) for different values of Q , as a functionof Hilbert space size |H| . Open symbols denote simulationswhere the Z Q symmetry is resolved, in which case the Hilbertspace size is the block size Q L − (we averaged data over all Q equivalent blocks). Filled symbols denote results for thefull Hilbert space of size Q L . The dashed lines represent thevalues of (cid:104) r (cid:105) GOE m obtained in Sec. III (taken from Table II),while (cid:104) r (cid:105) GOE is the numerical estimate for the GOE distri-bution taken from . The precision of our numerics allowsus to clearly distinguish the case m = 4 blocks with (cid:104) r (cid:105) GOE m =4 from the Poisson value (cid:104) r (cid:105) Poisson also represented in the plot.Simulations parameters are (cid:15) = 0 . J = 1, Γ = 0 . g = 0 . Q = 2, instead of the time-reversal and charge conjuga-tion breaking term in Eq. (61) which vanishes when Q = 2,we add a next-nearest neighbor interaction g (cid:80) j σ j σ j +2 inorder to break the mapping to a free-fermion model (we take g = 0 . Q eigen-states in the middle of the full spectrum, except for the smallersizes where ∼ Q eigenstates where considered. Results areaveraged over more than 4000 realizations of disorder, ex-cept for the largest size where 1000 realizations were used.For Q = 2 , , 4, we obtained results on chains of sizes up to L = 17 , , P ( r ), as obtained from simulations of chainsof sizes L = 16 , , Q = 2 , , P GOE m ( r ) obtained from the analyticalcomputations in Sec. III. B. Discrete symmetries in disorder distributions In this section, we consider the Heisenberg spin chainin the presence of a random external field h j : H Heisenberg = 12 L (cid:88) j =1 σ j · σ j +1 − L (cid:88) j =1 h j σ zj , (62)where σ x,y,z are the standard Pauli matrices, and we useperiodic boundary conditions. In general, this systemhosts a many-body localized phase at large enough dis-order: in particular, the model with box disorder has be-come the standard model of MBL in one dimension .When disorder is reduced, the system undergoes a tran-sition towards a thermal phase, a signature of which isan RMT-like spectral statistics. Since uncorrelated dis-order explicitly breaks all spatial symmetries, we expectthe spectral statistics in the thermal phase to be of GOEtype.However, in the specific case of binary disorder, i.e.taking on discrete values h j = ± h , visible deviationsfrom the GOE gap ratio distribution were observed inthe bulk of the thermal phase . The authors of ex-plained that this phenomenon was due to the peculiari-ties of discrete disorder distributions. Indeed, while ona finite-size system a typical disorder configuration willbreak all spatial symmetries, with discrete disorder dis-tributions such as the binary one, there is a non-zeroprobability that one or several of them are preserved,specially when considering periodic boundary conditions.For example, out of the 2 = 16 possible binary disorderconfigurations on L = 4 sites, 4 are reflection symmet-ric: (+ h, + h, + h, + h ), (+ h, − h, − h, + h ), . . . (actually, alldisorder configurations on L = 4 sites have a spatialsymmetry: reflection, translation, inversion – exchang-ing h ↔ − h , or a combination of them). Of course, when L is increased, the probability of drawing a spatially sym-metric configuration decreases exponentially fast; but forthe largest system sizes within reach of exact diagonaliza-tion techniques, the fraction of symmetric disorder con-figurations is still large enough to significantly alter thelevel statistics if disorder averaging is done “naively”,that is, by uniformly sampling over disorder. A possibleworkaround put forward in is to discard the spatiallysymmetric disorder configurations. If one insists on usingall samples, another possibility to is to explicitly resolvethe symmetry block structure of the Hamiltonian, when-ever the disorder configuration happens to be symmetric.This is cumbersome, especially given the number of pos-sible symmetries that must be taken into account.In that context, using the surmise for several GOEblocks proves useful. Indeed, the gap ratio in the ther-mal phase of the model can be written as a sum oversymmetry sectors: P ( r ) = (cid:88) m w m P GOE m ( r ) , (63)where w m is the weight of the symmetry sector of m . . . . . . r . . . P ( r ) . . . . . . δ P ( r ) FIG. 5. Probability distribution of the gap ratio from theHeisenberg model of Eq. (62) with L = 18 spins, with pe-riodic boundary conditions. Data is averaged over all re-alizations of the binary transverse field h j = ± / intotal, 3914 of which are nonequivalent up to symmetries),and over 150 eigenstates around infinite temperature energy E = ( E min + E max ) / 2. The red solid line represents theanalytical prediction Eq. (63): a linear combination of sur-mises P GOE m ( r ) obtained from the analytical computationsin Sec. III. For comparison, the green solid line shows thepredicted distribution when two blocks contributions are nottaken into account. Inset: Difference δP ( r ) between the nu-merical data and these surmises. blocks. In the thermodynamic limit L → ∞ , w → w m> decays exponentially fast to zero. Note thatfor samples with two blocks ( m = 2) the two blocks arealways of the same size, whereas for samples with morethan 2 blocks, m > 2, blocks are not necessarily of equalsize. While the expressions in the previous section allowus to compute the surmise for these non-homogeneoussamples, we can to a good degree of approximation ne-glect their contribution to the gap ratio distribution. In-deed, we find using simple combinatorics, that the totalweight (cid:80) m> w m coming from samples with more thantwo blocks is more than halved when L → L + 2, and for L = 18, it already represents less than 0 . 2% of the totalweight. We will therefore make the approximation that P ( r ) = w P GOE1 ( r ) + w P GOE2 ( r ).In Fig. 5, we show the gap ratio distribution for theHamiltonian Eq. (62) for L = 18. We find w =243936 / (cid:39) . w = 17640 / (cid:39) . 07. Thissystem size is of the order of what is achievable usingstate-of-the-art exact diagonalization techniques target-ing the middle of the energy spectrum . However, itis not large enough for w to be negligible comparedto w . Indeed, as shown in Fig. 5, incorporating the m = 2 contribution visibly improves the agreement withnumerical data. Note the clear difference at r = 0 be-tween Eq. (63) (for which P (0) (cid:54) = 0, as in the numer-ical simulations of Eq. (62)) and P GOE ( r ) which van-ishes at r = 0 due to level repulsion. Accordingly,the predicted average gap ratio using the m = 2 sur-mise (cid:104) r (cid:105) = w (cid:104) r (cid:105) + w (cid:104) r (cid:105) (cid:39) . 527 is closer to the numerically computed value (cid:104) r (cid:105) − (cid:104) r (cid:105) Heisenberg (cid:39) . m = 1: (cid:104) r (cid:105) m =1 − (cid:104) r (cid:105) Heisenberg (cid:39) . C. Floquet spin chain model We next consider both a static and a Floquet spin 1 / . Interacting many-body Floquet systems apriori exhibit no such interesting phases of matter, sincethe combination of interaction and driving is expected toheat up the system to an infinite-temperature, featurelessstate . However, it has been shown that the ad-dition of disorder, hindering energy propagation through-out the system via a MBL mechanism, can prevent heat-ing and give rise to new interacting Floquet phases, suchas discrete time crystals . Here, we study an interact-ing Floquet system, along with its static counterpart, forcomparison. We show in the following that the Floquetsystem exhibits an extra symmetry, that can be associ-ated to Floquet topological modes. In order to detect thesymmetry, we adjust the system parameters so as to bein the thermal phase of the model. Then, level statisticsis expected to follow RMT predictions, enabling us toemploy our surmises to detect the Floquet symmetry.We work with the following spin 1 / H x = L (cid:88) j =1 gσ xj H z = L (cid:88) j =1 Jσ zj σ zj +1 + L/ − (cid:88) k =0 h k +1 σ z k +1 , (64)again with periodic boundary conditions, which we com-bine to form a time-independent H static = H x + H z anda time-dependent model: H driven ( t ) = (cid:40) H z if 0 ≤ t mod τ < τ / H x if τ / ≤ t mod τ < τ . (65)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 isreplaced by quasi-energy , which is defined up to arbitraryshifts by 2 π/τ . More specifically, let us introduce the Flo-quet operator U F = exp( − iτ H x ) exp( − iτ H z ), which isthe evolution operator over one drive period. To the uni-tary Floquet operator we can associate a Floquet Hamil-tonian H F , defined as U F = exp( − iτ H F ), whose eigen-values ε α and associated eigenvectors are respectively thequasi-energies and the Floquet eigenstates, which hold4information about the dynamics and steady-state prop-erties of the system . In practice, when computing levelstatistics, we will therefore use the quasi-energies ε α ex-actly like energies in the static case.Going back to the system Eq. (64), remark that therandom longitudinal fields h j (cid:54) = 0 break both the Isingand the translation symmetries. Our model differs fromthe most commonly used one in that h j = 0 on evensites. This does not change the physics of the model,but can induce an extra symmetry in the driven case, aswe discuss below. Finally note that driving the systemdoes not break its time-reversal symmetry. We there-fore expect a GOE (respectively COE) level statistics inthe time-independent (respectively driven) case .Since the COE and GOE ensembles are asymptoticallydescribed by the same statistics, we will compare sim-ulations in the driven case to the corresponding GOEstatistics.We choose the parameter set g = Γ × . h k +1 =0 . 809 + 0 . × √ − Γ (cid:15) k , Γ = 0 . τ = π/ 4, wherethe (cid:15) k are uniformly distributed random number of zeromean and unit variance. This choice of parameters hasbeen shown to give good agreement with COE levelstatistics on the accessible system sizes, for the relatedmodel where the longitudinal field is also non-zero on theeven sites: h k (cid:54) = 0. This is indeed the case for the time-independent system, as can be seen in Fig. 6. However,the top panel of Fig. 6 shows that there is a dip in theaverage gap ratio (cid:104) r (cid:105) around J = 1 for the driven sys-tem. The numerical estimate for (cid:104) r (cid:105) at this special pointappears to coincide with the surmise value for two GOEblocks of equal size given in Table I. The level statis-tics (bottom panel of Fig. 6) is also compatible with thegap ratio statistics P GOE m =2 ( r ). Driving the otherwise fullyGOE system therefore appears to give rise to a new Z dynamical symmetry at the J = 1 point.We now rationalize the emergence of this symmetry.We find that the associated conserved operator is X = L/ − (cid:89) k =0 σ x k +1 . (66)Indeed, while this operator acts in general non-triviallyon the H z Hamiltonian, we have at the special point J =1, Xe iH z Xe − iH z = L (cid:89) j =1 iσ zj σ zj +1 = ( − L/ , (67)and thus X commutes with the Floquet operator, up to aglobal phase factor which can be absorbed in the defini-tion of U F . This indicates the existence of 2 COE blocksof identical sizes in the Floquet Hamiltonian, and henceexplains the agreement with the analytically-obtainedgap probability distribution P GOE2 ( r ). Note that setting h j = 0 on odd (or even) sites is necessary for this ex-tra symmetry to exist. Remark that when open bound-ary conditions are used, the above commutator Eq. (67) . 90 0 . 95 1 . 00 1 . 05 1 . J . . h r i h r i GOE h r i GOE2 h r i static h r i driven . . . . . . r . . . . P ( r ) P GOE1 ( r ) P static ( r ) P GOE2 ( r ) P driven ( r ) FIG. 6. Top: Average gap ratio as a function of J in thetime-independent and driven spin-1 / J = 1 for the time-independent and driven case (points), and comparison withthe surmise distribution for m = 1 and m = 2 GOE blocks(full lines). In the time-independent case, average is per-formed over 2000 disorder realizations (except at the J = 1point where 5000 realizations are used), and 100 eigenstatesaround infinite temperature energy E = ( E min + E max ) / 2, forthe system of size L = 14, while in the driven case, average isperformed over 2000 disorder realizations and all eigenstatesof the system of size L = 12. becomes proportional to σ z σ zL , a non-trivial boundaryterm. We can interpret this boundary term as creatinga pair of excitations . If the model were brought to theMBL regime (e.g. by increasing the strength of the dis-order term h j ), these excitations would become localizedat both ends of the chain, signalling the topological na-ture of the observed Z symmetry. However, in that casethe level statistics would become Poisson, and we wouldnot be able to detect the symmetry using our RMT ap-proach. Finally, we note that the argument carries overwhen we add a third contribution H y = (cid:80) L/ − k =0 h y σ y k +1 .It breaks the time-reversal symmetry, turning the 2-blockCOE structure of the J = 1 point into a 2-block CUEstructure.5 D. Anyonic Chain Our final application deals with chains of interactinganyons which are exotic particles interpolating betweenbosonic and fermionic statistics predicted to occur insome two-dimensional systems such as fractional quan-tum Hall states , and to offer exciting perspectives fortopological quantum computation . More precisely, wewill consider a disordered version of the “golden chain”model of Fibonacci anyons . As more technical back-ground is needed to introduce the model and its variousrepresentations, we first give a summary of our results.When periodic boundary conditions are imposed on thechain of anyons, there is a non-trivial topological sym-metry that decomposes the Hamiltonian into two blocksof unequal size. Resolving this symmetry is not easy,but it can in principle be done at the price of turningthe Hamiltonian into a dense matrix, rendering numeri-cal simulations on large systems difficult. Our results in-stead allow to use a representation simpler for numerics(with sparse, real symmetric matrices) which can nev-ertheless be confronted to RMT predictions, and henceprobe ergodic physics. This can be seen in Fig. 7 wherethe results for (cid:104) r (cid:105) (bottom curve) and P ( r ) allow to char-acterize the spectral statistics of the model with the twointerlaced sectors. At an extra numerical cost and withthe further requirement to study different representationsof the model, we can identify the states in each of thesetwo sectors, and they follow regular single-block GOEstatistics (two top curves in the top panel of Fig. 7). Inthe following, we present in detail the different repre-sentations of the model of disordered Fibonacci anyons,which allow us to draw these conclusions.The statistics of anyons are generally encoded in a setof fusion rules analogous to the composition rules for an-gular momenta, as well as transformation rules relatingthe different possible ways to fuse together three or moreparticles (the so-called “F-symbols”) . In the case of Fi-bonacci anyons there are only two types of particles, thetrivial particle, labeled by 1, and the Fibonacci anyon,labeled by τ , with fusion rule τ × τ = 1 + τ . This means,just as bringing together two spins 1/2 can be decom-posed into a spin 0 and a spin 1, that bringing togethertwo Fibonacci anyons can be decomposed as the directsum of a trivial particle, and another Fibonacci anyon.In addition we have 1 × × τ = τ . Followingthe seminal work , a Hamiltonian can be constructed byputting L particles of type τ along a chain and assigninga different energy for each possible outcome of the fusionof two nearest neighbors, H = − (cid:88) j J j Π j,j +1 , (68)where Π j,j +1 is the projector into the trivial particle oftwo τ particles located at sites j and j + 1. This isthe analog of the Heisenberg coupling for SU (2) spins1/2, which assigns a different energy to the fusion chan-nels of pairs of neighboring spins. The coupling con- stants J j are taken from a random distribution P ( J ) = (cid:15) − J − /(cid:15) θ ( J ) θ (1 − J ) ( (cid:15) ∈ [0 , ∞ ] characterizes the dis-order strength and J ∈ [0 , L anyons does not have a tensor product struc-ture: rather, it is convenient to describe it using a ba-sis of “fusion paths”. These correspond to sequences | x x . . . x L (cid:105) , where for each i , x i ∈ { , τ } , and x i +1 be-longs to the fusion τ × x i (this means in practice that notwo consecutive 1 are allowed in the sequence).In the case of periodic boundary conditions, x L +1 = x , the number of basis states |H| is related to the Fi-bonacci sequence, H ( L ) = F L − + F L +1 , where F L is the L th Fibonacci number with F = 0 and F = 1. It iswell known that the ratio of consecutive Fibonacci num-bers tends to lim i →∞ F i +1 /F i = φ with φ = √ thegolden mean, hence the name golden chain. A pedagog-ical introduction to the Hilbert space and Hamiltonianconstruction of the golden chain model can be found in .A practical representation for numerical simulations isto recast the above fusion paths in terms of sequencesof heights | h h . . . h L (cid:105) , where h i ∈ { , , , } and | h i − h i +1 | = 1, through the mapping 1 , τ → , i odd,1 , τ → , i even. This defines a “restricted solid-on-solid” (RSOS) model, namely the p = 4 case of the A p (also known as SU (2) p − ) family, where for generic p the heights are allowed to run between 1 and p .In this formulation the projectors Π j,j +1 of Eq. (68) canbe re-expressed in terms of operators e j , whose action isdefined as e j | h . . . h j − h j h j +1 . . . h L (cid:105) = δ h j − ,h j +1 (cid:88) h (cid:48) j (cid:113) sin( πh j p +1 ) sin( πh (cid:48) j p +1 )sin( πh j +1 p +1 ) | . . . h j − h (cid:48) j h j +1 . . . (cid:105) . (69)The operators e j form a representation of the Temperley-Lieb (TL) algebra , namely e j = √ Qe j , e j e j ± 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 oneindeed checks that Π j,j +1 = √ Q e j . Up to the irrelevant1 / √ Q proportionality factor, the Hamiltonian Eq. (68)is therefore re-expressed in the RSOS representation as H RSOS = − (cid:88) j J j e j . (70)A subtlety to keep in mind is that the RSOS formulationacts separately on two equivalent sectors, which corre-spond to putting even or odd heights on even sites re-spectively. For periodic boundary conditions, h L +1 = h ,each of these sectors has size F L +1 + F L − , and yields acopy of the original anyonic chain. The spectrum of theoriginal Hamiltonian Eq. (68) is therefore obtained by6restricting to a single of these sectors (which is what wedo in the following).The Hamiltonian Eq. (70) is real, a reason for whichthis representation is often used in numerics. We presentin Fig. 7 the results for the average gap ratio (cid:104) r (cid:105) and itsdistribution P ( r ) for Eq. (70), for different chain sizes(and thus Hilbert space sizes H ( L )) and weak disorder (cid:15) = 0 . 2, for states located in the middle of the spectrum.For this small value of disorder, we do expect a randommatrix theory behavior, but the value of (cid:104) r (cid:105) is clearlydifferent from the GOE statistics for a single block. Thesize of the Hilbert space hints at the existence of twoblocks of different sizes (denoted N and N in the fol-lowing). A first simple test is to compare the expecta-tion value (cid:104) r (cid:105) RSOS (cid:39) . 452 to the values obtained for twoGOE blocks of different sizes (Fig. 3 in Sec. III D). Thisleads to a possible value around α = N N + N ∈ [0 . , . N /N ∈ [0 . , . φ − = lim L →∞ F L − F L +1 (cid:39) . N = F L − and N = F L +1 . In the top panel of Fig. 7, we also rep-resent the predicted value (cid:104) r (cid:105) m =2 ,α =1 / (1+ φ ) = 0 . P ( r )(bottom panel of Fig. 7) which is an excellent match withthe one obtained from the surmise (see Sec. III) of twoGOE blocks with ratio φ − .We can in fact trace back this decomposition tothe existence of a “hidden” symmetry of a topologicalnature , namely an operator Y corresponding toan extra τ particle circling around the system, and whosematrix elements in the basis of fusion paths may be writ-ten as (cid:104) x (cid:48) . . . x (cid:48) L | Y | x . . . x L (cid:105) = L (cid:89) i =1 (cid:16) F x (cid:48) i +1 τx i τ (cid:17) x (cid:48) i x i − , (71)where the F-symbols (cid:16) F x (cid:48) i +1 τx i τ (cid:17) x (cid:48) i x i +1 can be found for in-stance in . The operator Y has two distinct eigenvalues (1 ±√ Y maps one onto the other. We can overcome thisdifficulty by computing the action of Y , which acts sep-arately in the odd and even sectors: this allows to definein each sector two orthogonal subspaces of dimensions F L +1 and F L − respectively, which precisely reproducesthe numerical observations made above .Now, it is important to remark that the action of Y ishighly non-local, and its matrix expression in the RSOSrepresentation is not sparse. Therefore while we knowin principle how to decompose the Hamiltonian into twoGOE blocks, it is not possible to our best knowledge to doso while keeping it sparse and real. One may ask whetherother representations of our model might help with this |H| . . . . . . h r i h r i GOE h r i GOE2 ,α =1 / (1+ φ ) Y = −√ Y = √ RSOS . . . . . . r . . . . . . . . P ( r ) P GOE2 ,α =1 / (1+ φ ) ( r )RSOS FIG. 7. RSOS model — Top: Average gap ratio (cid:104) r (cid:105) for theRSOS model Eq. (70), as a function of Hilbert space size.We distinguish between the results for the full spectrum, andthose in the Y = ±√ sectors. The dashed lines represent theresults for the (single-block) GOE distribution and from thesurmise computations in Sec. III for m = 2 GOE blocks withsize ratio φ − . Data in the Y = √ sector were obtained bycomparing energy spectrum in the loop representation (withnon-contractible loop weight 2 cos π/ 5) and RSOS representa-tion. Data in the Y = −√ were obtained by considering therest of the states in the RSOS representation. We focus onmid-spectrum eigenstates of the RSOS spectrum, obtaining ∼ 300 eigenstates for every disorder realization. Results areaveraged over between 300 and 1000 realizations of disorder ofstrength (cid:15) = 0 . 2. Bottom: Probability distribution of the gapratio P ( r ), as obtained from simulations of a RSOS chain ofsize L = 22. The solid line is the surmise P GOE m =2 ,α =1 / (1+ φ ) ( r )obtained from the analytical computations in Sec. III. problem. There are indeed other ways to represent theTL algebra, from which the spectrum of Eq. (68) can berecovered. Below we consider two such representations,the loop representation and the spin chain representation. Loop representation In the loop representation , theHilbert space is spanned by the configurations of non-crossing valence bonds between L vertical strands, andthe TL generator e i acts by contracting together thestrands at site i and i + 1. The composition rules ofthe TL algebra express the fact that lines can be contin-7uously deformed without crossing, and that closed loopscontribute a weight √ Q . From there, one can recoverthe eigenvalues of the anyon chain corresponding to eachsymmetry sector by assigning a special weight to non-contractible loops which close around the cylinder , re-spectively 2 cos π and 2 cos π , which is nothing but thecorresponding eigenvalue of Y . However, the loop modelcontains significantly more states than the anyonic chain,as the loops carry additional non-local information whichis absent in the RSOS representation. This brings severaldifficulties, the first being that the maximum size L thatcan be reached using exact diagonalization techniques islower, the second being that it is not obvious at all how toextract from the loop model spectrum the set of eigenval-ues which are present in the RSOS one . Furthermore,the loop representation leads to a non-Hermitian matrixrepresentation of the Hamiltonian, which also decreasesthe efficiency of simulations. Spin chain representation Another representation isin terms of a spin 1 / C ) ⊗ L ,on which the TL generators act as e i = − (cid:16) e i ϕL σ + i σ − i +1 + e − i ϕL σ − i σ + i +1 + cos γ σ zi σ zi +1 − − i sin γ σ zi − σ zi +1 ) (cid:19) . (72)Here the matrices σ x,y,zi act as Pauli matrices on the i th spin, and as identity elsewhere, and γ is defined by √ Q = 2 cos γ . The role of the twist parameter ϕ is analogto that of the weights of non-contractible loops in the geo-metrical representation. More precisely, the Hamiltonianbuilt out of Eq. (72) commutes with the global magne-tization S z = (cid:80) i σ zi , and the eigenvalues of the RSOSmodel are recovered in the S z = 0 sector upon setting ϕ = π and ϕ = π , respectively. As for the loop case,the S z = 0 sector however contains more states than theRSOS ones, leading to the difficulties mentioned above(see for other occurrences in related models). More-over this Hamiltonian is complex in the σ z basis, whichalso leads to a decreased numerical efficiency.We use simulations both in the loop and spin chain rep-resentations and checked on small systems (up to L = 18)that all states in the RSOS representation can indeed befound in the loop representation (using non-contractibleloop weight taking either 2 cos π or 2 cos π values) or thespin chain representation (using a twist taking either π/ π/ π are simpler (as allloops have the same weight) and we could reach largersystems. This allowed us to identify all states in the Y = √ sector for chains of size up to L = 24 (seeFig. 7).Besides allowing to identify this two-block structure in H RSOS chains with periodic boundary conditions, the ac-tual value of (cid:104) r (cid:105) and distribution P ( r ) for N /N = φ − will be useful as a marker of an ETH/ergodic phase whenincreasing the value of disorder. Indeed, it has been pro-posed that disordered SU(2) chains could lead to a new form of non-ergodic, critical, phase which behav-ior is different from a many-body localized phase. Thisputative new critical phase could be identified by the de-parture of spectral statistics from the references valuesdisplayed above. V. SUMMARY, RELATION TO PREVIOUSWORKS AND PERSPECTIVES In summary, we analyzed and computed the statisticsof the gap ratio r , an essential tool in diagnosing many-body quantum chaos, when the existence of symmetriesresults in a block structure of the matrix under consider-ation. The analytical results we obtain, based on an ex-tension of a seminal work of Rosenzweig and Porter , arevirtually indistinguishable from numerical simulations onlarge random matrices. While a closed form can onlybe obtained in limited cases, our formulation, based onEqs. (27),(28) and (8)), is compact and generic enough tobe implemented easily for all cases of interest. Throughseveral examples of applications, we showed the validityand usefulness of our results to identify or probe symme-tries in many-body quantum physics. In this final part ofthis manuscript, we relate our findings to previous works(including a re-interpretation of results available in theliterature) and provide leads for possible extensions. A. Relation to, and re-interpretation of previousresults We now relate our findings to others obtained in stud-ies of spectral statistics in various contexts. Our resultsallow to indirectly discover symmetries in a many-bodychaotic system, or to bypass them when they are toocomplex/costly to implement numerically. There havebeen several cases of unusual values of P ( r ) or (cid:104) r (cid:105) re-ported in previous literature which our work directly elu-cidates. For instance, it applies to the spectral statisticsof the Hamiltonian of the fractional quantum Hall effectwhen orbital inversion is not resolved in the numerics .Our analysis also explains the results obtained on the2d square lattice quantum Ising model in momentumsectors k = (0 , 0) and k = ( π, π ) where not all symme-tries were resolved. The value P ( r = 0) (cid:39) . Z symmetry there. Our analy-sis also accounts for the results in the one-dimensional t − t (cid:48) − V clean fermionic model of when the inversionsymmetry-breaking field is small. Another context whereour work is relevant is the bosonic SYK model with two-body interactions where the gap ratio distribution (seeFig. 6 in ) appears to be close to P GUE m =2 ( r ), suggesting atwo-block GUE structure (for instance due to a particle-hole symmetry), instead of an integrability signature asoriginally suggested in .The emergence of a symmetry not present in the origi-nal Hamiltonian may arise at critical points, with several8documented examples at quantum phase transitions (seee.g. ). While we are not aware of an example of anemerging symmetry for eigenstates deep in the spectrum,we can hypothesize that this could happen in the contextof transitions within or to a MBL phase. For instance,we note that the 2-block average value (cid:104) r (cid:105) GOE2 of the gapratio for the GOE/COE is extremely close to the onereported for a class of Floquet quantum circuits (withlocal qubit dimension 2) which behavior is between er-godic and MBL phases as well as for the FloquetMBL paramagnet / spin-glass transition in . While thismay be a pure numerical coincidence, it could specula-tively suggest an effective Z symmetry emerging at thesecritical points (we note however that recent works rather point towards an intermediate thermal phase in-stead of a direct transition between MBL phases). Ingeneral, transition points between an ergodic and a MBLphase systematically give to rise to values of (cid:104) r (cid:105) locatedbetween the single-block ergodic and Poisson values, butcloser to the latter (see e.g. ), with probability distri-bution showing a weak level repulsion. The method de-scribed in our work could allow to count the number ofeffective ergodic blocks or “bubbles” at or close to theMBL transition, albeit it would require a more detailedmodeling (as several combinations of blocks of differentsizes can give rise to the same expectation value (cid:104) r (cid:105) ). Ina different context, a similar computation could estimatethe number of partial barriers in the chaotic spectrum ofpressure modes in rapidly rotating stars, as studied in .In another direction, our analysis could be useful todiscover fracton models where the Hamiltonian de-composes in several different Krylov independent blocks(and this not necessarily based on an unresolved sym-metry), which necessarily implies a non-adherence to thesingle-block gap ratio statistics . B. Perspectives Our work can be extended in several directions. Ina straightforward way, it is possible to extend the anal-ysis to several blocks with different spectral statistics,for instance, the coexistence of GOE and GUE blocks inthe same spectrum. This applies to the quantum Hallwork of where different momentum sectors have differ-ent spectral statistics. Also, it is possible to see the effect of combining integrable and chaotic blocks, in the spiritof the work of on mixed phase spaces. This would applyto models with integrable “Krylov” subspaces co-existingwith ergodic blocks , or to the effective model of theMBL transition proposed in with one ergodic blockand random independent energies.Our approach is general enough that it should applymutatis mutandis to other RMT ensembles. Here we con-sidered the Wigner ensembles with quadratic potential inEq. (1). More generally, β -ensembles with different po-tentials can be treated in the same way. For instance, β -Laguerre ensembles, with logarithmic potential, are con-nected with Wishart matrices and are relevant to char-acterize entanglement spectra (for a review see e.g. ).Entanglement spectra can also exhibit block structureinherited from the symmetry of the underlying quantumstate from which they are formed. Other natural ex-tensions include replacing Eq. (3), which is our startingpoint, by the equivalent expression for matrices of largersizes: indeed obtains, from the exact expression of thejoint spacing distribution for 4 × P ( r ) which is more accurate than Eq. (7) by an orderof magnitude. Another possible direction is to study thenon-Hermitian situation with symmetries. Finally, itwould be interesting to analyze the case of weak symme-try breakings (with small matrix elements between dif-ferent blocks), using a perturbative approach to estimate P ( r ) in the same vein as the computation performed forthe level spacing distribution in . ACKNOWLEDGMENTS We thank B. Georgeot, L. Herviou for useful discus-sions as well as R. Vasseur for earlier insightful corre-spondence on the two-block structure of the periodicRSOS chain. This work benefited from the support of theproject THERMOLOC ANR-16-CE30-0023-02 and theproject COCOA ANR-17-CE30-0024-01 of the FrenchNational Research Agency (ANR), and by the FrenchProgramme Investissements d’Avenir under the pro-gram ANR-11-IDEX-0002-02, reference ANR-10-LABX-0037-NEXT. We acknowledge the use of HPC resourcesfrom CALMIP (grants 2018-P0677 and 2019-P0677) andGENCI (grant x2019050225). We use the librariesPETSc , SLEPc and Strumpack for themany-body computations presented in this manuscript. ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] D. J. Gross, The role of symmetry in fundamental physics,Proceedings of the National Academy of Sciences ,14256 (1996). E. P. Wigner, Characteristic vectors of bordered matrices with infinite dimensions, Annals of Mathematics , 548(1955). T. Guhr, A. MllerGroeling, and H. A. Weidenmller,Random-matrix theories in quantum physics: commonconcepts, Physics Reports , 189 (1998). M. L. Mehta, Random matrices (Elsevier, 2004). F. J. Dyson, Statistical theory of the energy levels of com-plex systems. i, Journal of Mathematical Physics , 140 (1962). M. Berry and M. Tabor, Level clustering in the regularspectrum, Proc. R. Soc. Lond. A , 375 (1977). O. Bohigas, M. J. Giannoni, and C. Schmit, Characteriza-tion of chaotic quantum spectra and universality of levelfluctuation laws, Phys. Rev. Lett. , 1 (1984). F. Evers and A. D. Mirlin, Anderson transitions, Rev.Mod. Phys. , 1355 (2008). A. Pal and D. A. Huse, Many-body localization phasetransition, Phys. Rev. B , 174411 (2010). F. Alet and N. Laflorencie, Many-body localization:An introduction and selected topics, Comptes RendusPhysique , 498 (2018). G. Montambaux, D. Poilblanc, J. Bellissard, and C. Sire,Quantum chaos in spin-fermion models, Phys. Rev. Lett. , 497 (1993). T. C. Hsu and J. C. Angl`es d’Auriac, Level repulsion inintegrable and almost-integrable quantum spin models,Phys. Rev. B , 14291 (1993). H. Bruus and J.-C. Angl`es d’Auriac, Energy level statis-tics of the two-dimensional hubbard model at low filling,Phys. Rev. B , 9142 (1997). V. Oganesyan and D. A. Huse, Localization of interactingfermions at high temperature, Phys. Rev. B , 155111(2007). B. Dietz and F. Haake, Taylor and pade analysis of thelevel spacing distributions of random-matrix ensembles,Zeitschrift f¨ur Physik B Condensed Matter , 153 (1990). F. Haake, Quantum Signatures of Chaos , 3rd ed., SpringerSeries in Synergetics (Springer-Verlag, Berlin Heidelberg,2010). J. M. G. G´omez, R. A. Molina, A. Rela˜no, and J. Re-tamosa, Misleading signatures of quantum chaos, Phys.Rev. E , 036209 (2002). Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux,Distribution of the ratio of consecutive level spacings inrandom matrix ensembles, Phys. Rev. Lett. , 084101(2013). Y. Y. Atas, E. Bogomolny, O. Giraud, P. Vivo, andE. Vivo, Joint probability densities of level spacing ratiosin random matrices, Journal of Physics A: Mathematicaland Theoretical , 355204 (2013). E. Cuevas, M. Feigel’man, L. Ioffe, and M. Mezard, Levelstatistics of disordered spin-1/2 systems and materialswith localized Cooper pairs, Nature Communications ,1128 (2012). D. J. Luitz, N. Laflorencie, and F. Alet, Many-body local-ization edge in the random-field heisenberg chain, Phys.Rev. B , 081103 (2015). V. Khemani, C. R. Laumann, and A. Chandran, Sig-natures of integrability in the dynamics of rydberg-blockaded chains, Phys. Rev. B , 161101 (2019). Y.-Z. You, A. W. W. Ludwig, and C. Xu, Sachdev-ye-kitaev model and thermalization on the boundary ofmany-body localized fermionic symmetry-protected topo-logical states, Phys. Rev. B , 115150 (2017). T. Kanazawa and T. Wettig, Complete random matrixclassification of SYK models with N = 0, 1 and 2 su-persymmetry, Journal of High Energy Physics , 50(2017). T. Li, J. Liu, Y. Xin, and Y. Zhou, Supersymmetric SYKmodel and random matrix theory, JHEP , 11 (2017). E. Iyoda, H. Katsura, and T. Sagawa, Effective dimension,level statistics, and integrability of sachdev-ye-kitaev-like models, Phys. Rev. D , 086020 (2018). F. Sun and J. Ye, Periodic table of the ordinary and su-persymmetric sachdev-ye-kitaev models, Phys. Rev. Lett. , 244101 (2020). P. Roushan, C. Neill, J. Tangpanitanon, V. M. Bastidas,A. Megrant, R. Barends, Y. Chen, Z. Chen, B. Chiaro,A. Dunsworth, A. Fowler, B. Foxen, M. Giustina, E. Jef-frey, J. Kelly, E. Lucero, J. Mutus, M. Neeley, C. Quin-tana, D. Sank, A. Vainsencher, J. Wenner, T. White,H. Neven, D. G. Angelakis, and J. Martinis, Spectroscopicsignatures of localization with interacting photons in su-perconducting qubits, Science , 1175 (2017). B. Evano, B. Georgeot, and F. Ligni`eres, Correlations inthe chaotic spectrum of pressure modes in rapidly rotatingstars, EPL (Europhysics Letters) , 49002 (2019). S. C. L. Srivastava, A. Lakshminarayan, S. Tomsovic,and A. Bcker, Ordered level spacing probability densities,Journal of Physics A: Mathematical and Theoretical ,025101 (2018). S. H. Tekur, U. T. Bhosale, and M. S. Santhanam, Higher-order spacing ratios in random matrix theory and complexquantum systems, Phys. Rev. B , 104305 (2018). U. T. Bhosale, S. H. Tekur, and M. S. Santhanam, Scal-ing in the eigenvalue fluctuations of correlation matrices,Phys. Rev. E , 052133 (2018). U. T. Bhosale, Superposition and higher-order spacing ra-tios in random matrix theory, arXiv:1905.02585 (2019). S. Harshini Tekur and M. S. Santhanam, Symmetry de-duction from spectral fluctuations in complex quantumsystems, arXiv:1808.08541 (2018). S. H. Tekur, S. Kumar, and M. S. Santhanam, Exact dis-tribution of spacing ratios for random and localized statesin quantum chaotic systems, Phys. Rev. E , 062212(2018). L. S´a, P. Ribeiro, and T. Prosen, Complex spacing ratios:A signature of dissipative quantum chaos, Phys. Rev. X , 021019 (2020). N. Rosenzweig and C. E. Porter, ”repulsion of energy lev-els” in complex atomic spectra, Phys. Rev. , 1698(1960). M. V. Berry and M. Robnik, Semiclassical level spac-ings when regular and chaotic orbits coexist, Journal ofPhysics A: Mathematical and General , 2413 (1984). A. J. Friedman, R. Vasseur, A. C. Potter, andS. A. Parameswaran, Localization-protected order in spinchains with non-abelian discrete symmetries, Phys. Rev.B , 064203 (2018). A. Prakash, S. Ganeshan, L. Fidkowski, and T.-C. Wei,Eigenstate phases with finite on-site non-abelian symme-try, Phys. Rev. B , 165136 (2017). F. Sun, Y. Yi-Xiang, J. Ye, and W.-M. Liu, Classifica-tion of the quantum chaos in colored sachdev-ye-kitaevmodels, Phys. Rev. D , 026009 (2020). E. P. Wigner, Statistical properties of real symmetric ma-trices with many dimensions (Princeton University, 1957). L. D’Alessio and M. Rigol, Long-time behavior of isolatedperiodically driven interacting lattice systems, Phys. Rev.X , 041048 (2014). The reader should not get confused by the fact that weinverted the notations between r and ˜ r with respect toRefs. . I. Dumitriu and A. Edelman, Matrix models for betaensembles, Journal of Mathematical Physics , 5830(2002). F. Pietracaprina, N. Mac´e, D. J. Luitz, and F. Alet, Shift-invert diagonalization of large many-body localizing spinchains, SciPost Physics , 045 (2018). P. Fendley, Parafermionic edge zero modes in z n -invariantspin chains, Journal of Statistical Mechanics: Theory andExperiment , P11020 (2012). P. Fendley, Free parafermions, Journal of Physics A:Mathematical and Theoretical , 075001 (2014). E. Fradkin and L. P. Kadanoff, Disorder variables andpara-fermions in two-dimensional statistical mechanics,Nuclear Physics B , 1 (1980). J. Alicea and P. Fendley, Topological phases withparafermions: theory and blueprints, Annual Review ofCondensed Matter Physics , 119 (2016). G. Albertini, B. M. McCoy, and J. H. H. Perk, Level cross-ing transitions and the massless phases of the superin-tegrable chiral potts chain, Physics Letters A , 204(1989). R. J. Baxter, The superintegrable chiral potts model,Physics Letters A , 185 (1988). E. Vernier, E. O’Brien, and P. Fendley, Onsager symme-tries in $u(1)$ -invariant clock models, Journal of Statis-tical Mechanics: Theory and Experiment , 043107(2019). J. Janarek, D. Delande, and J. Zakrzewski, Discrete dis-order models for many-body localization, Phys. Rev. B , 155133 (2018). A. Eckardt, Colloquium: Atomic quantum gases in pe-riodically driven optical lattices, Rev. Mod. Phys. ,011004 (2017). L. D’Alessio and M. Rigol, Long-time behavior of isolatedperiodically driven interacting lattice systems, Phys. Rev.X , 041048 (2014). A. Lazarides, A. Das, and R. Moessner, Equilibrium statesof generic quantum systems subject to periodic driving,Phys. Rev. E , 012110 (2014). P. Ponte, A. Chandran, Z. Papi, and D. A. Abanin, Peri-odically driven ergodic and many-body localized quantumsystems, Annals of Physics , 196 (2015). P. Ponte, Z. Papi´c, F. Huveneers, and D. A. Abanin,Many-body localization in periodically driven systems,Phys. Rev. Lett. , 140401 (2015). A. Lazarides, A. Das, and R. Moessner, Fate of many-body localization under periodic driving, Phys. Rev. Lett. , 030402 (2015). D. A. Abanin, W. De Roeck, and F. Huveneers, Theoryof many-body localization in periodically driven systems,Annals of Physics , 1 (2016). F. Harper, R. Roy, M. S. Rudner, and S. Sondhi, Topologyand broken symmetry in floquet systems, Annual Reviewof Condensed Matter Physics , 345 (2020). N. Y. Yao, A. C. Potter, I.-D. Potirniche, and A. Vish-wanath, Discrete time crystals: Rigidity, criticality, andrealizations, Phys. Rev. Lett. , 030401 (2017). V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi,Phase structure of driven quantum systems, Phys. Rev.Lett. , 250401 (2016). J. H. Shirley, Solution of the schr¨odinger equation witha hamiltonian periodic in time, Phys. Rev. , B979(1965). H. Kim, T. N. Ikeda, and D. A. Huse, Testing whether alleigenstates obey the eigenstate thermalization hypothesis,Phys. Rev. E , 052105 (2014). L. Zhang, H. Kim, and D. A. Huse, Thermalization of entanglement, Phys. Rev. E , 062128 (2015). L. Zhang, V. Khemani, and D. A. Huse, A floquet modelfor the many-body localization transition, Phys. Rev. B , 224202 (2016). T. L. M. Lezama, S. Bera, and J. H. Bardarson, Apparentslow dynamics in the ergodic phase of a driven many-bodylocalized system without extensive conserved quantities,Phys. Rev. B , 161106 (2019). W. Berdanier, M. Kolodrubetz, S. A. Parameswaran, andR. Vasseur, Floquet quantum criticality, Proceedings ofthe National Academy of Sciences , 9491 (2018). N. Read and E. Rezayi, Beyond paired quantum hallstates: Parafermions and incompressible states in the firstexcited landau level, Phys. Rev. B , 8084 (1999). H. Bartolomei, M. Kumar, R. Bisognin, A. Marguerite,J.-M. Berroir, E. Bocquillon, B. Plaais, A. Cavanna,Q. Dong, U. Gennser, Y. Jin, and G. F`eve, Fractionalstatistics in anyon collisions, Science , 173 (2020). C. Nayak, S. H. Simon, A. Stern, M. Freedman, andS. Das Sarma, Non-Abelian anyons and topological quan-tum computation, Reviews of Modern Physics , 1083(2008). A. Feiguin, S. Trebst, A. W. W. Ludwig, M. Troyer, A. Ki-taev, Z. Wang, and M. H. Freedman, Interacting anyonsin topological quantum liquids: The golden chain, Phys.Rev. Lett. , 160409 (2007). P. H. Bonderson, Non-Abelian Anyons and Interferom-etry , Ph.D. thesis, California Institute of Technology(2007). S. Trebst, M. Troyer, Z. Wang, and A. W. W. Ludwig, AShort Introduction to Fibonacci Anyon Models, Progressof Theoretical Physics Supplement , 384 (2008). V. Pasquier, Etiology of IRF models, Communications inMathematical Physics , 355 (1988). G. E. Andrews, R. J. Baxter, and P. J. Forrester,Eight-vertex sos model and generalized rogers-ramanujan-type identities, Journal of Statistical Physics , 193266(1984). H. N. V. Temperley and E. H. Lieb, Relations betweenthe percolation and colouring problem and other graph-theoretical problems associated with regular planar lat-tices: Some exact results for the percolation problem, Pro-ceedings of the Royal Society of London A: Mathematical,Physical and Engineering Sciences , 251 (1971). C. Gils, E. Ardonne, S. Trebst, D. A. Huse, A. W. W.Ludwig, M. Troyer, and Z. Wang, Anyonic quantum spinchains: Spin-1 generalizations and topological stability,Phys. Rev. B , 235120 (2013). M. Buican and A. Gromov, Anyonic chains, topologicaldefects, and conformal field theory, Communications inMathematical Physics , 1017 (2017). J. Belletˆete, A. M. Gainutdinov, J. L. Jacobsen, H. Saleur,and T. S. Tavares, Topological defects in periodic RSOSmodels and anyonic chains, arXiv:2003.11293 (2020). On a more algebraic note, the decomposition of the RSOSspace into F L +1 and F L − can be understood as a decom-position into irreducible representations of the periodicTemperley-Lieb algebra . More recently, it has been ob-served that the two orthogonal subspaces of the periodicRSOS model can be obtained from the open chain withfixed boundary conditions, through an operation called braid translation ? . H. Saleur and M. Bauer, On some relations between lo-cal height probabilities and conformal invariance, Nuclear Physics B , 591 (1989). E. Vernier, J. L. Jacobsen, and H. Saleur, Elaborating thephase diagram of spin-1 anyonic chains, SciPost Phys. ,004 (2017). V. Pasquier and H. Saleur, Common structures betweenfinite systems and conformal field theories through quan-tum groups, Nuclear Physics B , 523 (1990). H. Saleur, The antiferromagnetic potts model in two di-mensions: Berker-kadanoff phase, antiferromagnetic tran-sition, and the role of beraha numbers, Nuclear PhysicsB , 219 (1991). J. L. Jacobsen and H. Saleur, The antiferromagnetic tran-sition for the square-lattice Potts model, Nuclear PhysicsB , 207 (2006). E. Vernier, J. L. Jacobsen, and J. Salas, Q-colourings ofthe triangular lattice: Exact exponents and conformalfield theory, Journal of Physics A: Mathematical and The-oretical , 174004 (2016). R. Vasseur, A. C. Potter, and S. A. Parameswaran, Quan-tum criticality of hot random spin chains, Phys. Rev. Lett. , 217201 (2015). M. Fremling, C. Repellin, J.-M. St´ephan, N. Moran, J. K.Slingerland, and M. Haque, Dynamics and level statisticsof interacting fermions in the lowest landau level, NewJournal of Physics , 103036 (2018). R. Mondaini, K. R. Fratus, M. Srednicki, and M. Rigol,Eigenstate thermalization in the two-dimensional trans-verse field ising model, Phys. Rev. E , 032104 (2016). C. Cheng and R. Mondaini, Many-body delocalizationwith random vector potentials, Phys. Rev. A , 053610(2016). A. Nahum, P. Serna, J. T. Chalker, M. Ortu˜no, andA. M. Somoza, Emergent so(5) symmetry at the n´eelto valence-bond-solid transition, Phys. Rev. Lett. ,267203 (2015). C. Wang, A. Nahum, M. A. Metlitski, C. Xu, andT. Senthil, Deconfined quantum critical points: Symme-tries and dualities, Phys. Rev. X , 031051 (2017). H. Suwa, A. Sen, and A. W. Sandvik, Level spectroscopyin a two-dimensional quantum magnet: Linearly dispers-ing spinons at the deconfined quantum critical point,Phys. Rev. B , 144416 (2016). S. Gazit, F. F. Assaad, S. Sachdev, A. Vishwanath, andC. Wang, Confinement transition of z gauge theories cou-pled to massless fermions: Emergent quantum chromody-namics and so(5) symmetry, Proceedings of the NationalAcademy of Sciences , E6987 (2018). A. Chan, A. De Luca, and J. T. Chalker, Spectral statis-tics in spatially extended chaotic quantum many-bodysystems, Phys. Rev. Lett. , 060601 (2018). N. Mac, Quantum circuit at criticality, arXiv:1912.09489(2019). R. Sahay, F. Machado, B. Ye, C. R. Laumann, and N. Y.Yao, Emergent ergodicity at the transition between many-body localized phases, arXiv:2008.08585 (2020). S. Moudgalya, D. A. Huse, and V. Khemani, Perturbativeinstability towards delocalization at phase transitions be-tween mbl phases, arXiv:2008.09113 (2020). M. Pretko, X. Chen, and Y. You, Fracton phases ofmatter, International Journal of Modern Physics A ,2030003 (2020). P. Sala, T. Rakovszky, R. Verresen, M. Knap, and F. Poll-mann, Ergodicity breaking arising from hilbert space frag-mentation in dipole-conserving hamiltonians, Phys. Rev. X , 011047 (2020). V. Khemani, M. Hermele, and R. Nandkishore, Localiza-tion from hilbert space shattering: From theory to phys-ical realizations, Phys. Rev. B , 174204 (2020). S. Moudgalya, A. Prem, R. Nandkishore, N. Regnault,and B. A. Bernevig, Thermalization and its absencewithin krylov subspaces of a constrained hamiltonian,arXiv:1910.14048 (2019). X. Wei, R. Mondaini, and G. Xianlong, Characteriza-tion of many-body mobility edges with random matrices,arXiv:2001.04105 (2020). S. N. Majumdar, Extreme eigenvalues of wishart matri-ces: application to entangled bipartite system, in TheOxford Handbook of Random Matrix Theory , edited byP. D. F. G. Akemann, J. Baik (Oxford University Press,2011) Chap. 37. D. M. Leitner, Real symmetric random matrix ensemblesof hamiltonians with partial symmetry breaking, Phys.Rev. E , 2536 (1993). S. Balay, W. D. Gropp, L. C. McInnes, and B. F. Smith,Efficient management of parallelism in object oriented nu-merical software libraries, in Modern Software Tools inScientific Computing , edited by E. Arge, A. M. Bruaset,and H. P. Langtangen (Birkh¨auser Press, 1997) pp. 163–202. S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune,K. Buschelman, L. Dalcin, V. Eijkhout, W. D. Gropp,D. Kaushik, M. G. Knepley, L. C. McInnes, K. Rupp,B. F. Smith, S. Zampini, H. Zhang, and H. Zhang, PETScUsers Manual , Tech. Rep. ANL-95/11 - Revision 3.8 (Ar-gonne National Laboratory, 2017). V. Hernandez, J. E. Roman, and V. Vidal, SLEPc: Ascalable and flexible toolkit for the solution of eigenvalueproblems, ACM Trans. Math. Software , 351 (2005). J. E. Roman, C. Campos, E. Romero, and A. Tomas, SLEPc Users Manual , Tech. Rep. DSIC-II/24/02 - Revi-sion 3.8 (D. Sistemes Inform`atics i Computaci´o, Univer-sitat Polit`ecnica de Val`encia, 2017). P. Ghysels, X. Li, F. Rouet, S. Williams, and A. Napov,An efficient multicore implementation of a novel HSS-structured multifrontal solver using randomized sampling,SIAM J. Sci. Comput. , S358 (2016). P. Ghysels, X. S. Li, C. Gorman, and F. H. Rouet, A ro-bust parallel preconditioner for indefinite systems usinghierarchical matrices and randomized sampling, in2017IEEE International Parallel and Distributed ProcessingSymposium (IPDPS)