The universal spectral form factor for many-body localization

We theoretically study universal correlations present in the spectrum of many-body-localized systems. We obtain an exact analytical expression for the spectral form factor of Poisson spectra and show that it agrees well with numerical results on two models exhibiting a many-body-localization: a disordered quantum spin chain and a phenomenological l-bit model based on the existence of local integrals of motion. We find that the functional form of the Poisson spectral form factor is distinct from but complementary to the universal expectation of quantum chaotic systems obtained from random matrix theory.

The understanding of how thermal equilibrium may, or may not emerge in isolated many-body quantum systems remains a central question in quantum statistical mechanics. Thermal systems which are said to exhibit quantum chaos satisfy the eigenstate thermalization hypothesis (ETH) [1,2] whose subsystems equilibriate under their own dynamics. In addition to being highly entangled, i.e. "volume law", the eigenstates of these systems exhibit long range repulsions that are captured by random matrix theory and produce universal features in their spectral form factor (SFF) such as the linear ramp [3][4][5][6] (as shown in fig. 1). In the presence of strong quenched randomness or quasiperiodicity, quantum systems can become many body localized (MBL) [7][8][9][10] where ETH is violated. In contrast to chaotic systems, MBL is characterized by short-range "area law " entanglement and an absence of level repulsion. Recent experiments on ultra-cold atomic gases [11][12][13], trapped ions [14], superconducting qubits [15,16] and nuclear spins [17] have provided evidence for the existence of the MBL phase.
Instabilities to MBL have been argued to arise in high dimensions [18] and in the presence of certain symmetries [19]. More recently however, the very existence of the MBL phase itself has been challenged based on a finite size scaling analysis of the linear ramp of the SFF on approach to the MBL transition from the chaotic side [20]. A rebuke to this work has subsequently been presented [21] pointing out the intricacies involved in finite sized calculations and conclusions drawn from them, while further studies have highlighted the difficulty in studying the MBL transition in finite size numerics [22,23]. Recently, the authors of Ref. [20] pointed out that their claim of the absence of MBL is due to their choice of scaling function, which instead should follow a Kosteritz-Thouless "like" scaling form as they demonstrate in Ref. [24] consistent with recent theories of the MBL transition [25][26][27]. Irrespective of the question of validity of the finite-size numerics in the vicinity of the MBL transition, the question of how to characterize the MBL phase using the SFF alone, is undoubtedly worthy of further examination. If the MBL phase indeed  The adjacent gap ratio, r defined in eq. (3) versus W . The approximate critical disorder, where the data at different system sizes cross is given by Wc ≈ 7.3 has also been marked. For W Wc in the MBL phase the level statistics are Poisson r = 2 ln(2) − 1 ≈ 0.39 [29]. The dashed black line is the well known GOE expectation from random matrix theory known to describe the thermal phase, whereas the solid black line, that matches the numerical data in the MBL phase well (over the range of W ≥ 10 [30]), is obtained in this work in eq. (8).
exists, it is conceivable that its SFF has its own universal functional form and features to which any putative system exhibiting MBL should be compared with. However, apart from a few hints [28], the existence of such a universal form and its features has been lacking thus far.
In this letter, we investigate the universal spectral correlations in MBL systems. We show that, like quantumchaotic systems, the SFF for MBL systems indeed has a universal form due to spectra consisting of Poisson levels for which we derive an exact analytical expression (plotted in solid line in fig. 1). We determine the validity of our result by comparing it with numerical simulations of a phenomenological l-bit model [31,32] as well as a microscopic disordered many-body Hamiltonian. In both cases, by focusing on states in the middle of the many body spectra where the many-body density of states is nearly flat, we find excellent agreement between our exact expression and the numerical results. Our results provide further support for the existence of the MBL phase in one dimension.
Models for many-body localization: To make a detailed comparison with the properties of the MBL phase, we consider two different models.
The first is the Hamiltonian of a quantum spin-chain with quenched disorder, which is defined as where S α are spin operators that can be written in terms of Pauli matrices as S α = 1 2 σ α and the random couplings w i are drawn from a uniform distribution [−W, W ]. Variants of this model have been previously studied in Refs [20,33,34], which are known to have a thermal phase at weak disorder and an MBL phase at strong disorder. Following Ref 20, we set J 1 = J 2 = 1.0 and ∆ = 0.55.
Deep in the MBL phase, any local Hamiltonian such as eq. (1) can be described by a complete set of emergent local integrals of motion [31,32]. As a result, there should exist a finite depth circuit of local unitary operators U that can recast H into a diagonal form, U HU † = H lbit : (2) where κ z i are the so called l-bit Pauli operators with localized support on the Hilbert space near site i, whose eigenvalues represent the locally conserved quantities and the magnitudes of J m i1...im fall off exponentially with distance. The second model we consider is a truncated version of the above phenomenological l-bit model eq. (2) where we retain only up to 10 spin nearest neighbor interactions with all couplings drawn from a uniform distribution J (1...10) ∈ [−1, 1].
Characterizing spectral correlations of quantum systems: A popular diagnostic used to distinguish MBL and chaotic systems via their spectral correlations is the adjacent gap ratio (r) [29]. This is defined in terms of successive gaps δ i = E i+1 − E i of an ordered energy spectrum {E i } as follows: For chaotic systems, the value of r (where . . . denotes averaging over samples and energy) can be computed from an appropriate random matrix ensemble. For example, for systems with time-reversal symmetry, the Gaussian orthogonal ensemble (GOE) gives r ≈ 0.53, while for MBL systems, r = 2 ln(2) − 1 ≈ 0.39 [29]. As shown in the inset of fig. 1, by tracking r , we can see that the Hamiltonian of eq. (1) supports a thermal phase for small W and an MBL phase for large W with the critical disorder strength somewhere near W c ≈ 7.3 where the curves for different system sizes cross, consistent with previous work [20]. The adjacent gap ratio captures the repulsion of neighboring levels, and thus only probes local spectral correlations. It does not measure long-range spectral correlations which has important and useful information. A more powerful diagnostic is the spectral form factor (SFF) [3], which is the main tool of analysis in this letter and is defined in terms of eigenvalues {E i } as follows where N is the total number of eigenvalues in consideration. We can also define the connected SFF , The information about long-range correlations is contained in the form of K(τ ) interpolating the fixed early and late τ values of N 2 and N respectively (0 and N for K c (τ ) respectively). For chaotic systems, just like r , the SFF can also be computed from an appropriate random matrix ensemble. For instance, as seen in fig. 1, the SFF for the Hamiltonian eq. (1) for small disorder W exhibits a clear ramp and matches that of the GOE ensemble whose approximate expression (plotted in dotted line) is known [3,5,35,36]: where J 1 is the Bessel function of the first kind. As we increase disorder W , as shown in fig. 1, the SFF qualitatively changes as the model passes through the MBL transition with the disappearance of the ramp. Deep in the MBL phase (i.e. where r ≈ 0.39), the SFF again takes on a new universal form K P (τ ) (plotted as a solid black line in fig. 1). The expression for K P (τ ) as well as the connected version, K P c (τ ) we obtain are [30] We will show in the following that eqs. (8) and (9) correspond to energy levels drawn from a Poisson process. Note that both expressions K GOE (τ ) and K P (τ ) as well as the data in fig. 1 are normalized to set the mean level spacing to unity. Contrasting features between K GOE and K P can be seen at intermediate τ values, in the regime where the SFF is non-trivial (this occurs in the range 1/D τ 1/µ where D = µN is the many-body bandwidth of the chosen levels and µ is the mean level spacing [3,5,6,35,36]). For K GOE , this corresponds to the 'ramp region' governed by the connected piece, K GOE c . This has an N dependence only by an overall scaling and does not affect the slope of the ramp on a log-scale. On the other hand, as expected, both K P and the connected version K P c lack the ramp. They also depend quite sensitively on N as can be seen in fig. 2

(dashed lines).
Random levels versus Poisson levels: While the single particle spectrum of the non-interacting Anderson insulator [37] can simply be described by a set of random levels, the same is not true about the corresponding many body spectrum. The latter is not a single but a weighted sum of random numbers. In the presence of interactions, however, how this perspective changes is unclear, and it is not obvious how to describe the spectra of MBL systems.
To get a hint, we recall that deep in the middle of the MBL spectrum where the many-body density of states is nearly flat, it is has been observed that the spectral gaps δ i = E i+1 − E i of N levels follow an exponential independent and identical distribution (i.i.d) [38] where µ is the mean level spacing. This distribution of gaps reproduces the observed value of the adjacent gap ratio r ≈ 0.39. The computation of the SFF eqs. (4) and (5) however requires not just the gaps but the distribution of the full spectrum {E i } which is not uniquely determined from that of the gaps δ i alone.
If we regard the gap distribution eq. (10) to be a necessary condition, there are in fact at least two distinct possibilities for the full spectrum {E i }: (1) The spectrum of a random diagonal matrix (that we refer to as random levels). (2) Numbers drawn from a Poisson process (that we refer to as Poisson levels). The SFF is calculated exactly for these two cases below and shown to be fundamentally different. In what follows, with no loss of generality, we will assume that the energy eigenvalues are bounded from below and non-negative E ≥ 0.
First we consider the case (1) i.e. N ×N random diagonal matrices with entries all drawn from some i.i.d, p(E). The spectrum, which is basically the sorted N random numbers on the diagonal, has gaps whose distribution is known to converge to eq. (10) for large N [39]. The SFF for this case, which we denote as K R (τ ) can be easily

computed [30] as
where, e iτ E is computed from the i.i.d p(E). For example, if the numbers are drawn from a uniform distribution The SFF for this case is plotted in fig. 2 in solid lines where we have set the mean level spacing to unity by fixing D = N/2.
Next, we consider the case (2) i.e. eigenvalues drawn from a Poisson process known to describe the distribution of radio-active decay events [40], The above spectrum is directly obtained by summing up from the gaps {δ j } of eq. (10). Once again, we can derive the SFF in a straight forward manner [30], which yields eqs. (8) and (9) stated previously (the dependence on mean level spacing µ can be restored by replacing τ → µτ ) For clarity, we note that the term 'Poisson' is sometimes also used to refer to the distribution of gaps eq. (10) rather than the levels eq. (13). To avoid confusion, we will reserve the term only for the levels and refer to the gap distribution in eq. (10) as exponential.
A few comments are in order: (i) The SFF for the two cases are distinct from each other as seen in fig. 2. At intermediate τ values, K P (τ ) 'bounces off' and K P c (τ ) 'hovers around' 2N which is twice the asymptotic value of N to which both 'roll off' at large τ . In contrast, as can be seen in eq. (12), The shape of the SFF, i.e the number of oscillations in K P c and bounces in K P , is set by the number of levels N , which is in strong contrast to the thermal regime where N enters as a prefactor. (iii) At this stage, it is not obvious which (if either) of the two cases is relevant to MBL. We answer this in the next section by numerical analysis and find that the second case is the correct one.
Comparison with numerical calculations: We now numerically check if either of the above possibilities are applicable to the spectrum of physical systems in the MBL phase. We consider the two models defined in eq. (1) (disordered quantum spin chain) and eq. (2) (phenomenological l-bit model). Both models possess a global U (1) spin rotation symmetry which allows us to focus on half-filling i.e. the total S z = 0 sector. We will perform our analysis by shifting the N chosen eigenvalues by the smallest one so as to make them non-negative. For ease of comparison with the analytical results as well as across system sizes, after averaging over disorder samples, we re-scale τ by the mean level spacing µ effectively setting the mean level spacing to unity. We stress that N is chosen independently of the system size L for the computation of K(τ ). Depending on system size, our analysis is performed using disorder samples ranging from 5, 000 to 50, 000 [30].
It is a well known challenge to compare exact random matrix theory predictions with numerics on microscopic models due to the latter possessing a non-uniform density of states that produces non-universal corrections [5].
To deal with this, we focus deep in the middle of the spectrum where the non-universal corrections from curvature in the many-body density of states (DOS) is minimized. At finite sizes however, deviations are bound to exist, e.g. even if the middle of the DOS is flat the tails on the edges of the many-body density of states will strongly affect the results at short τ (that probe large energy differences) in a non-universal way. Nonetheless, these corrections to our analytic results are expected to reduce by either increasing the system size L and thus (2) (above) and the microscopic Hamiltonian in eq. (1) deep in the MBL phase with W = 25 (below) for various system sizes (L) and eigenvalues (N ) that are compared with the analytical curves K P (τ ) of eq. (8). We find excellent agreement in the limit of vanishing ζ = N/NL, where NL is the total Hilbert space size. For increasing ζ we see deviations appear that we attribute to non-universal structure in the many-body density of states that is not captured by our theory [30].
the total Hilbert space size N L = L L/2 (due to working in the total S z = 0 sector) or choosing a smaller set of eigenvalues N . In the thermodynamic limit (L → ∞) when the parameter ζ ≡ N/N L vanishes for any finite N , we expect any deviations to completely vanish and the analytical results to match exactly [30].
We start with the l-bit model of eq. (2). Since it is already diagonal, the eigenvalues are generated easily and as a result, we are able to reach relatively large system sizes. As seen in fig. 3 (top panel), the numerical SFF, K(τ ) matches the analytical one for Poisson levels, K P (τ ) of eq. (8) (dotted lines) very well with negligible deviations even for values of N as large as 1000.
We now turn to the microscopic Hamiltonian eq. (1) and focus deep in the MBL phase at W = 25, where r ≈ 0.39 is nicely Poisson at the accessible L. Here, we are relatively limited in the system sizes that we can reach and the presence of complex microscopic details further impacts the finite sized numerical results more severely than in the case of the idealized l-bit model. Nevertheless, as seen in fig. 3 (lower panel), for small values of N (20,40), the numerical SFF matches the analytical eq. (8) (dotted lines) very well. For larger values of N ∼ 80, deviations due to DOS curvature start becoming visible, which most notably occur at short τ values that are most strongly affected by the tails on the edges of the manybody density of states. The root mean square error of our exact result and the numerics systematically decreases as ζ → 0 [30]. Although we have only presented the analysis for W = 25, we find that the results remain virtually unchanged for a wide range of disorder strengths, with the universal SFF holding for W ≥ 10 [30]. This strongly supports the notion that MBL is a robust phase in disordered one-dimensional many-body systems.
Conclusion: In this letter, we have derived an exact universal expression for the spectral form factor of Poisson levels. We have shown that this describes the shape of the SFF in the many body localized phase through a detailed comparison with numerical results on two separate models that describe many-body localized systems.
Note added: During the final stages of drafting our paper, we became aware of a recent mathematical physics paper [41] which also briefly discusses the spectral form factor for spectra with uncorrelated spacings in a distinct context.

Supplementary Material
In this supplemental material we provide details on the calculations of the spectral form factor (SFF) for random and Poisson levels that we compare with a numerical evaluation. The many-body density of states is presented, as well as the system size dependence of the SFF both computed using the model Hamiltonians in the main text. The comparison of the universal SFF and data in the MBL phase for a range of disorder strengths is presented, showing remarkable agreement across the entire phase. Last, we analyse the root mean square error of our numerical results with the exact SFF prediction.

S1. SFF for random levels
We provide details of the derivation of the expressions for the spectral form factors (SFF) of random levels presented in the main text. Consider an N × N matrix with random entries along the diagonal drawn from some independent and identical distribution (i.i.d) p(E). The eigenvalues {E i } of such a matrix is simply the diagonal entries sorted. As can be easily seen from the formulas for the SFF, The ordering of {E i } is irrelevant to compute the SFF. We need the joint two-point distribution P (E n , n; E m , m), i.e. probability that the m th eigenvalue is E m and the n th eigenvalue is E n . Since we can ignore ordering, this is simply Let us start with the disconnected SFF for random levels, K R (τ ) where we have used e iτ Em e −iτ En P = e iτ Em p e −iτ En p , We can also get the connected SFF easily which gives us the expressions in the main text. Observe that K c (τ ) in the above expression strictly satisfies the inequality As an example, if the levels are drawn from a uniform distribution [0, 2D], the expressions for SFF are as follows (S11) Figure S1 shows a comparison of the analytic result for the SFF for random levels and with it, the numerically computed SFF for N random numbers drawn from a uniform distribution, E m ∈ [0, N ] (width selected to ensure unit mean level spacing) and its excellent matching with the above analytical expression. . The connected (left) and disconnected (right) SFF for random levels with unit mean level spacing. The analytic result is the black line for N = 500 levels, which we compare with the numerically computed SFF using levels drawn from a uniform distribution that is averaged over 50,000 realizations.

S2. SFF for Poisson numbers
We now provide details of the derivation of the expressions for the spectral form factors (SFF) of Poisson levels presented in the main text. Consider levels from a Poisson process {E i } generated by summing gaps {δ i } with an exponential i.i.d To compute the SFF, again, we need the joint two-point distribution P (E n , n; E m , m) i.e. the probability that the m th eigenvalue is E m and the n th eigenvalue is E n . Assuming with no loss of generality m > n, this can be obtained as follows where p(E k , k) is the well known Poisson distribution where x > 0 and n spans over non-negative integers 0, 1, 2, . . .. Let us sketch the derivation of the above equation where (S20) In the last step, the integral is performed by closing the contour in the upper-half complex plane to enclose the n th order pole at k = i/µ and using Cauchy's integral formula.
Let us now proceed to compute the SFF where where we have used ∞ 0 dE n p(E n , n) = 1. Substituting the above expression in eq. (S21) gives The summand of the above series depends only on the differences k = m−n. As a result, we can change the summation over m and n in terms of k as follows The above series can be summed up easily using standard results from arithmetico-geometric series S n = ab + (a + d)br + (a + 2d)br 2 + . . . + (a + (n − 1)d)br n−1 to get the final answer We now proceed to calculate the connected SFF where Plugging this into eq. (S28) and some simplification, we get the final answer (S31) Figure S2 shows the SFF computed numerically for N Poisson levels and the excellent agreement with the above analytical result. We now focus on the disorder dependence of the SFF in the model of the disordered spin-chain that we considered in the main text. We list it here again for completeness As described in the main text, by keeping track of the adjacent gap ratio r that is defined in terms successive gaps we can determine that the model above has a thermal phase for low disorder (W < W c ) and an MBL phase at high disorder (W > W c ). The critical disorder is roughly around W c ≈ 7.3 where the curves for different system sizes appear to cross as shown in fig. S3. In the main text, we focused on data that is deep in the MBL regime (i.e. where r ≈ 0.39 is nicely Poisson), for a specific disorder strength W = 25 where we showed that the SFF converges well with the analytical result. We now test the validity of this result across the MBL phase. Figure S4 shows that for a wide range of disorder strengths W ≥ 10, the SFF is universal and remarkably convergent with each other and, importantly, with the analytical result. This further strengthens the claim that MBL is indeed a robust phase.  As mentioned in the main text, at finite sizes, deviations from the exact result due to a non-uniform many-body density of states (DOS) is inevitable. We minimize this by picking out N eigenvalues deep in the middle of the spectrum where the DOS is closest to flat. To examine this, we compute the many-body DOS defined as for N eigenvalues i = E i /L. As shown in fig. S5, we observe that for N corresponding to increasingly smaller fractions of the total Hilbert space dimension in the S z = 0 sector, N L = L L/2 , the DOS appears increasingly flat for both the disordered spin chain model H(W ) of eq. (S32) in the MBL phase (W=25) as well as the truncated l-bit model eq. (S35) (we list here again the l-bit model for completeness).
Despite removing corrections from curvature in the DOS in the middle of the many-body spectrum, our results at small τ values are dominated by large energy differences, which probe the edges of the many-body spectrum. As shown in fig. S5, due to finize size there are tails on the edges of the many-body DOS that sharpen for increasing L and decreasing N . We attribute the deviations in our results at short τ values with the numerics to these band edge effects which contribute non-universal corrections that are not included in the analytic result. The DOS curvature and thus the deviation from the analytical result are expected to reduce by increasing system size (L and thus the size of the Hilbert space L L/2 ) and reducing N . This is clearest in the connected SFF as shown in fig. S6. Consistent with our expectations from the many-body DOS, we find the largest area of deviation occurs at short τ values, where as larger τ values (but well before the plateau "time") that are sampling energy differences from the flat DOS, our analytic results match the numerics nicely.
To quantify the deviation, we consider the quantity ∆ defined as the root mean square difference of the numerical connected SFF from the analytical value at various values of τ as defined below.
τ i is chosen from the interval τ i ∈ (10 −5 , 15.0) (with mean-level spacing set to 1) where the SFF is non-trivial and N τ is the total number of τ points. As seen in fig. S7, we see ∆ decrease as expected with increasing L as well as decresing N .