On intermediate statistics across many-body localization transition

The level statistics in the transition between delocalized and localized {phases of} many body interacting systems is {considered}. We recall the joint probability distribution for eigenvalues resulting from the statistical mechanics for energy level dynamics as introduced by Pechukas and Yukawa. The resulting single parameter analytic distribution is probed numerically {via Monte Carlo method}. The resulting higher order spacing ratios are compared with data coming from different {quantum many body systems}. It is found that this Pechukas-Yukawa distribution compares favorably with {$\beta$--Gaussian ensemble -- a single parameter model of level statistics proposed recently in the context of disordered many-body systems.} {Moreover, the Pechukas-Yukawa distribution is also} only slightly inferior to the two-parameter $\beta$-h ansatz shown {earlier} to reproduce {level statistics of} physical systems remarkably well.


I. INTRODUCTION
Fritz Haake made fundamental contributions first to quantum optics (see e.g. [1,2]) then to the fast developing in the eighties of the last century area of Quantum Chaos not only with his seminal monograph on the subject [3] but also with many original works from introducing a celebrated kicked top model [4] to providing a link between Gutzwiller's periodic orbit theory [5] and random matrix statistics [6][7][8].
Those late works provided a highlight of Haake's fascination of the link between spectral properties of physical models and random matrix theory. Of particular interest for him has always been level dynamics (see Chapter.6 of [3]). The dependence of the energy levels of a Hamiltonian H(λ) on some parameter λ may be viewed as the motion of interacting fictitious particles (levels) with "time" λ following the original formulation of Pechukas [9] and Yukawa [10]. For large matrices of size N × N , in N → ∞ limit, statistical mechanics is applicable to those fictitious particles. Haake devoted particular interest to periodically driven (Floquet) systems also in this context, in particular since a transition between level clustering for integrable systems to level repulsion (for quantally chaotic system) may be viewed as a relaxation toward equilibrium [11].
For this transition, following the original Wigner 2 × 2 matrices approach, Lenz and Haake found an interpolating spacing distribution [12] which compared well with random matrix simulations adding a significant contribution to the topic which originated with the early work of Rosenzweig and Porter [13]. The Rosenzweig-Porter distribution was an important twist on random matrices. As it is well known [3,14,15] for generalized time reversal invariant systems (the case we shall solely concentrate on) the Gaussian orthogonal ensemble of random matri- * jakub.zakrzewski@uj.edu.pl ces (GOE) may be constructed by considering matrices with independent Gaussian distributed entries with the variance of the diagonal elements being twice the variance of the off-diagonal elements. GOE ensemble faithfully represent statistical properties of quantally chaotic system as conjectured by Bohigas, Giannoni, and Schmidt [16]. The Rosenzweig-Porter ensemble interpolates between GOE and Poisson cases using a single parameter χ = σ 2 /N , where σ 2 is the variance of off-diagonal elements and N the matrix size [13]. Lenz and Haake solution [12] gives the nearest spacing distribution for such a 2 × 2 matrices ensemble. This is by no means a unique solution for the transition between GOE and Poisson level statistics. The different models that were proposed over the years are reviewed in Section II. The renewed interest in the problem comes from a realization that the transition between spectral properties of ergodic many-body interacting systems and the so called many-body localized phase [17][18][19] may be described in the statistical sense by an ensemble interpolating between GOE and Poissonian case as in the case of chaotic systems or single particle Anderson localization. In Section III we review the statistical approach of Yukawa [10] following the excellent textbook of H.-J. Stockman [15]. Such a Pechukas-Yukawa statistical mechanics model provides another, to the best of our knowledge not tested yet, proposition for the interpolating ensemble. This approach has certain aesthetic advantages over other propositions opening also possibilities of further studies going beyond the eigenvalue statistics. A critical comparison of this ensemble predictions with numerical data obtained for disordered Heisenberg chain and further examples for other systems are the content of Section IV. We conclude giving future perspectives of the presented approach in Section VI. arXiv:2108.11654v2 [cond-mat.dis-nn] 12 Nov 2021

II. THE INTERPOLATING ENSEMBLES AND FORMULAE
The early approaches to level statistics often concentrated, as the work of Lenz and Haake [12], on the nearest neighbor level spacing distribution. The distribution of normalized spacings s i = e i+1 −e i where e i are the Hamiltonian matrix eigenvalues renormalized (in the process called the level unfolding [3,15]) in such a way that their mean density is unity. Some proposed an ad hoc expressions as a popular Brody distribution [20] which surprisingly well fitted low resolution experimental data -see e.g. [15]. In contrast, Berry and Robnik [21] proposed a distribution that was based on a clear physical assumption of the separation between "chaotic" wavefunctions faithful to GOE and those localized on the regular parts of the phase space. This simple foundation resulted in the distribution which was shown to work well for mixed phase space situations in the deep semiclassical limit [22]. Importantly, long range correlations of energy levels can also be addressed within this model [23].
Over the years, the most common approach to construct an ensemble interpolating between GOE and uncorrelated Poisson level statistics has been to modify the GOE distribution resigning from its beautiful property of being invariant under orthogonal transformations. The original Rozenzweig-Porter proposition [13] is just one of them. In the original version a single parameter λ = σ 2 /N defined the ensemble of matrices of rank N where σ 2 was the variance of off-diagonal matrix elements of the Hamiltonian matrix. For λ → 0 one recovers the Poisson ensemble while for λ = 1 the GOE ensemble is reproduced provided that the variance of the diagonal elements is equal to 2. Recently, modifications of the Rosenzweig-Porter ensemble allowing for a broad, lognormal distribution of the off-diagonal matrix elements, have been shown to be valid for various types of manybody and Anderson localization transitions [24][25][26].
A second proposition goes back to Seligman and coworkers [27] who postulated that the variance of offdiagonal elements a ij should scale as exp −(i − j) 2 /σ 2 . For σ → 0 one recovers the Poisson case while for σ → ∞ the matrix belongs to GOE.
Another well-known approach is that of Guhr [28] who used supersymmetric techniques to express the twolevel correlation function in the Poisson-GOE ensemble in terms of a double integral. In the meantime Casati and coworkers [29,30] as well as Fyodorov and Mirlin [31] considered banded Gaussian random matrices as a useful tool in describing the transition, the corresponding parameter was y = b 2 /N with b being the matrix bandwidth and N its rank.
For a long time, the nearest neighbor spacings were treated as essential entities to capture the GOE-to-Poisson transition in the context of quantum chaos. An another relevant parameter is the number-variance, defined as Σ 2 (L) = n(L) 2 − n(L) 2 , where n(L) is the number of (unfolded) eigenvalues in an interval of length L. It was demonstrated that for a large system size L, the number variance scales as χL, where χ varies from 0 to 1 during the Possion-to-GOE transitions [32][33][34]. However, the unfolding process, implemented to compute e i and the number variance, is based on a separation of the density of states ρ(E) into smooth and fluctuating parts. The distinction between these two parts of ρ(E) is to some extent arbitrary and may lead to dubious outcomes see, for instance, the differences between local unfolding and Gaussian broadening [35]. Over the past years, there is a constant effort to tackle such problems from the mathematical and computational viewpoints [36,37].
To simplify the problem, Huse and Oganesyan [38] introduced an important measure known as gap ratio, defined as r n = min[δ n , δ n−1 ]/max[δ n , δ n−1 ], where δ n = E n+1 − E n is the energy gap between the consecutive energy levels. The dimensionless gap-ratio is independent of the local density-of-states, ρ(E), and the problem of unfolding is avoided. Employing an exact diagonalization, they also established that a spin system with time reversal invariance undergoes a phase transition from ergodic (with GOE level statistics) to the many-body localised (MBL) phase (signified by the uncorrelated Possonian statistics of energy levels) as the magnitude of disorder increases (see also [39]). In a further endeavor, Atas et.al. [40] proposed a Wigner-like surmise for the gap ratio statistics. While the level statistics across the ergodic-to-MBL transition was still studied with standard level spacing distributions [41], the gap ratio became gradually a more popular tool [42][43][44].
The quest for a proper level statistics model to describe ergodic-to-MBL transition speeded up following the work of Serbyn and Moore [45] who argued that the flow of level statistics between the GOE and Poisson limits occurs in two stages: (1) A power-law type interaction between the eigenvalues [46] of range that decreases with increase of disorder strength and (2) an appearance of semi-poissonian level statistics between the energy levels as an offshoot of local interactions between the energy levels (see also [32,47]).
Among the latest developments a so called β-Gaussian(β − G) model was introduced [48] as a model to describe level statistics across ergodic-to-MBL transition. The β − G model is dependent on a single parameter β that characterizes the pairwise interaction of the energy levels taking all the energy pairs in consideration. A real value of β ∈ [0, 1] allows one to describe the whole crossover between GOE and Poisson level statistics. For a comparison of the performance of different models see [49]. In our recent work [50] , we have proposed a modified 2-parameter model, namely, the β − h model, where the interaction between the energy levels is limited to a h neighboring energy levels. This 2-parameter model was shown to be superior to its β − G counterpart since it reproduced not only the nearest neighbor spacing ratio but also the higher order spacing ratios defined by [51,52]: where E i stands for the i-th energy level. In the present work we compare predictions of these recent models with the so called Pechukas-Yukawa (P − Y ) level-dynamics approach and verify its underlying efficiency to reproduce the level statistics in terms of the higher order spacing ratios (1). In the next section we remind the derivation of the joint-probability-distribution of the P − Y model (which may be found, e.g. in [15]) to make the paper self contained while later we compare the statistics generated with different models in the transition between delocalized and MBL cases.

III. STATISTICAL APPROACH TO LEVEL DYNAMICS
Consider a Hamiltonian dependent on parameter λ: H(λ) = H 0 + λV for arbitrary H 0 and V and the fate of its eigenvalues as λ is varied. Differentiating the eigenvalue equation with respect to λ (above x(λ) a is the eigenvalue corresponding to eigenvector |a(λ) ) and taking appropriate scalar products one immediately gets Defining p a ≡ẋ a and looking at its λ derivative yieldṡ where the second equality involves an additional definition f ab = V ab (x a − x b ). Most interestingly, differentiating f ab over λ, the set of differential equations closes: The set of equations (3)-(5) is known as the Pechukas-Yukawa equations. Identifying λ as a fictitious time, x a 's and p a 's can be interpreted as positions and momenta of fictitious particles moving under the force decaying as 1/x 3 with the distance (assuming f ab 's are constant). This highly nonlinear set of of equations is neverthless integrable (for a discussion see [3,15]). Clearly, for λ large the eigenvalues of H = H 0 + λV will be V dominated and the variation of eigenvalues (and eigenvectors) trivializes. For that reason Haake [3] introduces a slightly different λ dependence (equivalent for small λ): while we shall use the form adopted from [15,53,54] H = H 0 cos(λ) + V sin(λ). Then the equations above become slightly modified withẋ with f ab = a|Ḣ|b (x a − x b ). Thus the additional harmonic force binds the fictitious particles preventing their escape to infinity. Pechukas-Yukawa equations corresponding to integrable system call for an appropriate statistical mechanics description which should be realized within the generalized Giggs ensemble [55]. This requires an identification of the complete set of nontrivial constants of motion being in convolution. Such a possible approach is a song of the future. Here we rather follow Yukawa approach and consider the simplest integrals of motion. One of them will be the total energy of the system of interacting classical particles and the other Q = 1 2 T r(F 2 ) called a total angular momentum, for a justification see [15]. Let us note that those are the only two second order constants of the motion.
With this two constants the phase-space density is given by the Gibbs ensemble as which in the presence of the harmonic binding potential reads explicitly Next, we can compute the joint probability distribution (JPD) of eigenvalues [15] by integrating out the variables p n and f nm from the phase space density obtaining where ν = 1, 2, 4 corresponds to three possible ensembles of Dyson. All three cases appear due to different properties of (integrated over) f nm variables. We shall consider the generalized time reversal invariant case of ν = 1 only. From the expression above, one can easily reach the Possonian distribution by putting γ/β 1. The distribution becomes proportional to while in the opposite limit, 0 < γ/β << 1, the distribution mimics the Gaussian ensemble.
To obtain an ensemble that interpolates between GOE and Poisson level statistics we fix β = ν = 1 and denote γ/β = 10 p . The distribution (11) becomes where p = log 10 γ β is the single free parameter of the distribution. The first term in the (14) signifies the pairwise interaction between the particles and the exponential term provides the harmonic binding of the eigenvalues. This single parameter distribution, resulting from Pechukas-Yukawa statistical mechanics for time reversal invariant case, will be tested against other existing propositions as well as numerical data from interacting manybody models.

IV. COMPARISON WITH NUMERICAL DATA
For the simulation purpose, the eigenvalues are drawn by sampling (14) with Metropolis-Hastings algorithm [56,57] for p ∈ [−10, +15] with an initial sample size N = 500 distributed over a linear chain. This moderate system size is chosen for convenience and speed. We consider about 100 eigenvalues in the middle of the spectrum. Since we consider high order gap ratios too, this sample size may be insufficient for r n with n large. We consider the issue of the system size in the next Section. The obtained eigenvalue distributions interpolate between a faithful semi-circle histogram of the densityof-states in the GOE limit (when p → −∞) and a Gaussian distribution for p large and positive. In the following we extract the higher order spacing ratios from (14) and fit different physical models results with that distribution. For comparison we use the single parameter β-Gaussian model [48] as well as the two-parameter β − h model shown to yield accurate representation of many body data [50]. As a first example we take a disordered Heisenberg spin-1 2 chain, a paradigmatic model for MBL studies [18,44,45,49,58,59]. The Hamiltonian of the system readŝ where the first term represents the exchange interaction between the neighboring 1/2-spins with exchange coupling J (normalized to unity). The last term constitutes the on-site disorder potential at site i with disorder h i drawn from uniform random distribution, It is well established that this system undergoes a transition from ergodic to MBL phase [18,44,45,49,58] for sufficiently large W . The eigenvalues of the Heisenberg chain are obtained for the length L = 18 with periodic boundary conditions (PBC) via the shift-and-invert method. In the due course, we use 2000 disorder realizations for each value of W . The gap ratios are calculated from approximately 500 eigenvalues from the middle of the spectrum. In Fig. 1, we present a comparison between the obtained gap-ratio distributions P (r n ) and those resulting from the closest predictions coming from the P − Y model (14). Different n values allow us to explore correlations on different distances between eigenvalues. Let us stress that we use a single p value and the resulting P − Y distribution to fit all P (r n ) for n ∈ [1, 10] by minimizing the cumulative error between the gap-ratios of Heisenberg chain and the gap-ratios predicted by P − Y model. As may be appreciated in Fig. 1, the P − Y distribution fits the data of the spin model remarkably well, in particular on the delocalized side of the transition (smaller W values). Only close to the fully localized case some small deviations are seen. To see how competitive P − Y model is we compare its predictions with those given by a single parameter β − G approach, in which sampling of the distribution is much easier [48] as well as a two parameter family of β − h distributions [50]. For a most direct comparison with data in Fig. 1 we present similar figures for β −G in Fig. 2 and for β − h ensemble in Fig. 3. Even a casual inspection of Fig. 2 reveals that β − G ensemble fitted for r 1 (where it performs remarkably well) does not represent properly higher gap ratios. Clearly this ensemble does not reproducer short-range spacings and long-range correlations equally well. On the other hand β − h ensemble performs at least as well as the P-Y distribution considered.  For a global, quantitative comparison, we construct a measure of the fit quality, i.e., the average deviation where r n H is the numerical average for the Heisenberg chain while r n M is the corresponding value coming from a fit to a given model, be it P −Y , β −h or β −G version.
The obtained values of are plotted in Fig. 4(a), for β−G model β is fitted to reproduce as closely as possible r 1 H . Close to the delocalized GOE-like regime the resulting for that ensemble is comparable to both P − Y and β −h values. Closer to the transition and on the localized  side one may clearly observe that β − G model performs poorly when compared with P − Y and β − h models. The latter is only slightly superior but it involves fitting of two parameters instead of one. Fig. 4(b) presents the average gap ratios r n for n ∈ [1,40] or rather, for a better visualization, a difference ∆ r n = r n − r P n where the latter corresponds to an analytic prediction for the mean gap ratio for a Poissonian ensemble [50] (note that normalization factor is missing there). On a first glance it is clear, that the β − G model is inaccurate beyond r 1 as it severely overestimates ∆ r n , which means that the long-range correlations between eigenvalues of β − G model are much stronger than the correlations of eigenvalues of the disordered Heisenberg spin chain . Both β − h and P − Y reproduce higher order spacing ratios much more accurately. The single parameter P − Y model is only slightly inferior to the β − h model.

B. Quasi-periodic Heisenberg spin chain
Consider the same Heisenberg chain with the Hamiltonian (15) but now h i are not random but taken as h i = (W/2) cos(2πχi + α), where i is the site index. Different realizations of the disorder correspond to different choices of α drawn from random uniform distribution on [0, 2π) interval. Contrary to the previous case, this quasi-periodic disorder is fully correlated. Still, for sufficiently large W (dependent on the value of χ [60,61]) the system undergoes a transition to MBL phase. Such a quasiperiodic disorder was implemented in experiments of the Munich group [62,63] by placing an additional weak standing wave on top of the primary one forming the optical lattice. Then χ is the ratio of wavevectors of both laser wavelengths. In our study we take χ = ( √ 5 − 1)/2 and diagonalize the Hamiltonian for L = 16 and open boundary conditions. Like in the preceding case, we implement the shift-invert method to collect eigenvalues from 2000 disorder realizations and the higher order spacing ratios are computed form 500 eigenvalues taken from the middle of the spectrum.
The whole analysis is performed similarly to the random disorder case. We concentrate on the localized side of the crossover (for the system size studied) where the differences between the considered distributions show more clearly. Fig. 5 shows the resulting fits for P − Y ensemble showing that indeed it works very well also for this case. For comparison results obtained for β − G and β − h ensembles are shown in Fig. 6 and Fig. 7, respectively. The conclusions are very similar to the previous case (although limited to smaller interval of W taken). Firstly, the single parameter P − Y model outperforms β-Gaussian proposition considerably, especially for higher order gap ratios. Similarly to the uniform random disorder case, β-Gaussian model works reasonably well only very close to the GOE limit (data corresponding to W = 1.5 in Fig. 5. Secondly, P − Y distribution is clearly less effective than the phenomenological, two parameter β − h approach in the whole interval of W values. But, let us stress, this difference is quite small.
The global error comparison presented in Fig. 8 additionally confirms the above conclusions.
It is worthwhile to comment more on the comparison of random and quasiperiodic cases. The crossover to MBL present in both models is quite well reproduced by P − Y or β − h statistical models despite the fact that the systems with random and quasiperiodic disorder behave quite differently on a microscopic level. This was observed in [64] inspecting the entanglement entropy statistics and further analysed, on the level of level statistics, in our earlier work [49]. We have shown there that one may distinguish the MBL transition in random and quasiperiodic disorder by examining the so called intersample variances which are the sample-to-sample variances of gap ratios. Those show a peak in the crossover region for a random disorder case. This fact was attributed to the presence of rare Griffiths-type regions in systems with random disorder. The corresponding peak was absent in the quasiperiodic data supporting rare regions interpretation. One may pose a question: how inter sample variances of gap ratio behave in the statistical models such as P − Y or β − h models? We have verified that inter-sample variations of gap ratio practically do not change in the GOE-Poisson transition for the P − Y model by a direct evaluation. Our explanation of this discrepancy between behavior of a system with random disorder at the MBL transition and P − Y model is the following. The rare regions affect variances of the former, for which the number of independent random variables determining the Hamiltonian is the same as the system size, L, very small compared to the Hilbert space dimension. For the latter, the number of random variables in simulated matrices is much larger and scales as a square of the matrix size making the sample-to-sample variations very small across the whole GOE-Poisson crossover.
C. Bose-Hubbard model As the last example we test P − Y model on a different disordered system namely the Bose-Hubbard model with the Hamiltonian [64][65][66][67] where a † k (a k ) are bosonic creation and annhilation operators at site k, J is the tunneling and U is the onsite interaction strength. We assume further on J = U = 1 while the chemical potential µ i is assumed to be uniformly randomly distributed within an interval [−W, W ]. As the considered previously spin chain a transition from delocalized phase to MBL is observed for such a disor-dered Bose-Hubbard model with increasing W [65]. We consider a chain of L = 8 sites at unit filling and eigenvalues are computed via exact-diagonalization method [49]. As previously, the gap ratios are determined from the eigenvalues collected from the middle of the spectrum for no less than 500 disorder realizations. The distribution of gap ratio P (r n ) obtained are fitted with the P − Y model as shown in Fig. 9 for n = 1, 2, 3, 4, 5, 8. One may observe a rather outstanding agreement between the data and the statistical P − Y model except for slight deviations close to full MBL (Poisson) cases.
The corresponding comparison of numerical data for the disordered Bose-Hubbard model with β −G and β −h models are shown in Fig. 10 and Fig. 11, respectively. Similarly to previously considered spin chains , P − Y model is clearly superior over the β − G model in reproducing the data while being slightly outperformed by the two parameters β − h model, as further summarised in Fig. 12. The results presented for all three cases discussed show convincingly that the single parameter Pechukas-Yukawa distribution (14) reproduces remarkably well distributions of higher order spacing ratios in the transition between GOE and Poisson regime. Let us recall the fact that we have consistently used data for N = 500 ensem-ble. It is clear, from the form of the distribution (14) that in the Poisson limit the density of states becomes a Gaussian with a unit variance independent of the system size. In the other, GOE limit the density of states follows a semicircle law in (− √ 2N , √ 2N ) interval being strongly N dependent. In effect for any given set of data, the value of the fitted parameter p in (14) depends on N chosen.
This estetic drawback may be cured by rewriting (14) as where in the denominator the factor 10 p is replaced by 10 Y /∆ 2 , where ∆ is the mean spacing of {x n }. This defines a new single parameter Y characterizing the distribution, while ∆ = ∆(N, Y ) can be easily obtained from generated sequence of x n . We have verified that such an approach makes the obtained distributions N independent as visualised in Fig. 13 which presents the average gap ratios r i for i = 1, 3 for N = 500 and N = 2000.

VI. CONCLUSIONS
The aim of this work was to investigate how the distributions obtained from JPD of eigenvalues resulting from standard Pechukas-Yukawa statistical approach compare with currently discussed models of level statistics interpolating between GOE and integrable cases. The nice feature of P − Y model, as given by (14), is that a single parameter, denoted by p allows for an interpolation in the whole interval between GOE and Poisson limits. The distribution studied results in statistics closely following the numerical data and in this respect is clearly superior to other single parameter model such as β-Gaussian proposition [48]. The latter has a clear advantage that probing it does not require the Metropolis algorithm, so probing it is relative simple. We also made a comparison with the so called β −h model [50] which quite accurately represents numerical data. The difficulty of applying the latter are comparable. It is quite rewarding to see that a single parameter P − Y model yields results only slightly inferior to two-parameters β − h approach.
It is worth stressing that while we concentrated on eigenvalue statistics, the Pechukas-Yukawa JPD contains also information on distribution of matrix elements f nm of a specific operator F = [H,Ḣ]. The latter is manifestly time independent -thus matrix elements f nm carry information on eigenstates dependence on the parameter. Indeed Nakamura and Lakshmanan [68] rephrased the dynamics in terms of eigenvalues and components of eigenvectors. Thus JPD integrated over eigenvalues may yield prediction on eigenvector statistics in the transition regime (or alternatively on matrix elements of a generic operator).
Last but not least let us mention that one could aim at more "proper" analysis of integrable model given by (3)-(5) with the appropriate generalized Gibbs ensemble [55] provided a proper identification of all independent constants of the motion is made. A more modest approach will just include some of them going beyond the presented here quadratic terms approach.
The presented study may be generalized to cases with broken generalized time-reversal invariance as well as to an analogous approach possible for time-periodic Floquet systems.