Distributions of consecutive level spacings of Gaussian unitary ensemble and their ratio: ab initio derivation

In recent studies of many-body localization in nonintegrable quantum systems, the distribution of the ratio of two consecutive energy level spacings, $r_n=(E_{n+1}-E_n)/(E_{n}-E_{n-1})$ or $\tilde{r}_n=\min(r_n,r_n^{-1})$, has been used as a measure to quantify the chaoticity, alternative to the more conventional distribution of the level spacings, $s_n=\bar{\rho}(E_n)(E_{n+1}-E_n)$, as the former makes unnecessary the unfolding required for the latter. Based on our previous work on the Tracy-Widom approach to the Janossy densities, we present analytic expressions for the joint probability distribution of two consecutive eigenvalue spacings and the distribution of their ratio for the Gaussian unitary ensemble (GUE) of random Hermitian $N\times N$ matrices at $N\to \infty$, in terms of a system of differential equations. As a showcase of the efficacy of our results for characterizing an approach to quantum chaoticity, we contrast them to arguably the most ideal of all quantum-chaotic spectra: the zeroes of the Riemann $\zeta$ function on the critical line at increasing heights.

1. Introduction.Many-body localization that prohibits thermal equilibration of the wave functions has been a core agenda of research in the field of quantum many-body systems, including disordered and interacting fermions on a chain [1]; spin chains with transverse field, periodical kicks, or disorder [2][3][4][5]; and the Sachdev-Ye-Kitaev model and its deformed, coupled, or sparsed variant [6][7][8][9], to name but a few.In all of the above studies, the "gap ratio distribution," i.e. the distribution P r (r) of the ratio of consecutive energy level spacings r n = (E n+1 − E n )/(E n − E n−1 ) or rn = min(r n , r −1 n ), initiated by Ref. [1], has been utilized as a criterion of (non)ergodicity.This trend is obviously due to its computational advantage that makes unnecessary the unfolding by the smoothed density of states ρ(E), over the conventional distribution P (s) of level spacings s n = ρ(E n )(E n+1 − E n ).More recently, the use of the gap ratio distribution for characterizing spectral transitions extended its reach beyond energy level statistics of many-body or chaotic Hamiltonians, namely towards quantum entanglement spectra of reduced density matrices in neural network states [10] and in quantum circuits [11,12].
The random matrix theory of the gap ratio distribution was introduced in Ref. [2] and has since been quoted in practically every work in these fields including Refs.[3][4][5][6][7][8][9][10][11][12].There, two types of approximate expressions for P r (r) were presented: one from the Wigner-like surmise (substituting the large-N limit of N × N matrices by 3 × 3) for all three Dyson symmetry classes, P r (r) ≈ C β (r + r 2 ) β /(1 + r + r 2 ) 1+3β/2 (β = 1, 2, 4), and another from the quadrature discretization of the resolvent operator K I (I − K I ) −1 where K I denotes the integral operator of convolution with the sine kernel over an interval I for the unitary class β = 2.Although the latter approximation is known to converge to the exact value quickly as the quadrature order is increased [13], the analytical expression for P r (r) is still missing.Moreover, the author cannot help feeling frustrated to find that every article that quotes Ref. [2] almost always refers only to the Wigner-surmised form, and calls that crude, uncontrolled approximation the "outcome from the random matrix theory".In view of these, this article aims to determine analytically the joint distribution of consecutive eigenvalue spacings P c (a, b) and the distribution of their ratio P r (r) for the unitary class, based on our recent work [14] which provided a generic prescription for determining the Jánossy density for any integrable kernel as a solution to the Tracy-Widom (TW) system of partial differential equations (PDEs) [15].
This article is composed of the following parts: In Sect. 2 we list several facts on the Jánossy density J 1 (0; I) for the sine kernel, i.e. the conditional probability that an interval I in the spectral bulk of the GUE contains no eigenvalue except for the one at a designated locus.In Sect. 3 we follow the prescription of TW and show that the Jánossy density and associated distributions of consecutive eigenvalue spacings and their ratio are analytically determined as a solution to a system of ordinary differential equations (ODEs).As a showcase of the efficacy of our results for characterizing an approach to quantum chaoticity without using ρ(E), in Sect. 4 we contrast them to the distribution of zeroes of the Riemann ζ function on the critical line at increasing heights.A program for generating the Jánossy density is attached as Supplementary material.
2. Jánossy density for the sine kernel.The sine kernel governs the determinantal point process of unfolded * eigenvalues or eigenphases {x i } of random Hermitian or unitary matrices of infinite rank N → ∞, in the spectral bulk [16,17].
From the very defining property of the determinantal point process that the p-point correlation function R p (x 1 , . . ., x p ) = E[ p i=1 δ(x − x i )] is given by det[K(x i , x j )] p i,j=1 , it follows that the conditional p-point correlation function in which one of the eigenvalues is preconditioned at x = t, also takes a determinantal form det[ K(x i , x j )] p i,j=1 governed by another kernel [14], In the case of the sine kernel (1), its translational invariance allows for setting t = 0 without loss of generality, so that [18], The transformation of kernels (3) from K to K is associated with a meromorphic SL(2, R) gauge transformation [14], on the two-component functions that comprise respective kernels in the right-hand sides of Eqs. ( 1) and ( 4).Accordingly, as stated in Theorem in Ref. [14], the gauge-transformed section [φ(x), ψ(x)] T inherits from the original section [π −1/2 sin(x), π −1/2 cos(x)] T the covariant-constancy condition for a meromorphic sl(2, R) connection, i.e. a pair of linear differential equations, which guarantees applicability of the TW method.In the present case of spherical Bessel functions (4), the polynomials comprising the connection (Lax operator) Subsequently we shall use their nonzero coefficients, We take an interval I = [a 1 , a 2 ] with a 1 < 0 < a 2 so that the ordered triple (a 1 , 0, a 2 ) will serve as three consecutive eigenvalues, and denote by K I the integration operator acting on the Hilbert space of square-integrable functions L 2 (I) with convolution kernel K(x, y), Then by Gaudin and Mehta's theorem [16], the Jánossy density J 1 (0; [a 1 , a 2 ]) [19], i.e. the conditional probability that the interval I = [a 1 , a 2 ] contains no eigenvalue except for the one preconditioned at x = 0, is expressed as the Fredholm deterninant of K I : Note that the Jánossy density for a symmetric interval I = [−t, t] was previously expressed in terms of a Painlevé V transcendent, i.e. a special solution to an ODE in t [20].Our task in this article is to extend their result to a generic interval and express J 1 (0; [a 1 , a 2 ]) and an associated distribution P r (r) in terms of a system of PDEs in a 1 and a 2 .3. TW system.TW [15] established a systematic method of computing the Fredholm determinant of an integrable integral kernel whose component functions satisfy the condition 3/9 in Eq. ( 6).The quantities that appear in the TW system [j, k = 1 or 2], are all treated as functions of the left and the endpoints (a 1 , a 2 ) of the interval I.
Expanding the definitions in Eqs. ( 10) and ( 11) in a 1 and a 2 ≪ 1, the boundary conditions for these quantities read: up to terms of O(a 1 , a 2 ) 6 .Substituting the coefficients (8) into the TW system of PDEs (Eqs.(2.25), (2.26), (2.31), (2.32), (2.12)-(2.18),(1.7a) of Ref. [15]), it takes the following form [below the pair of indices (j, k) assumes either (1, 2) or (2, 1)]: where U = 1 + u − w and V = 2v.The "stiffness" of the second and the third equations of Eqs. ( 13) at a j = 0 is only superficial, because q j , p j , and R jk are of O(a j ) or higher orders 4/9 Fig. 1 The Jánossy density J 1 (0; [a 1 , a 2 ]) for the sine kernel (left), and the joint distribution P c (a 1 , a 2 ) of two consecutive eigenvalue spacings (right).For visual clarity, each coordinate a j is rescaled from the text by 1/π, so that the mean eigenvalue spacing is unity.
It is quite plausible that the ODEs ( 15) can be expressible as a Hamiltonian system, and J 1 (0; [a 1 , a 2 ]) be regarded as a τ -function of an integrable hierarchy [21].However, these are not immediately evident to the author, and will be discussed in a separate publication.As a cross-check of our formulation, we confirmed that the values of J 1 (0; [a 1 , a 2 ]) obtained by the above prescription (Fig. 1, left) are identical, within the accuracy of numerical evaluation, to those from the Nyström-type quadrature approximation of the Fredholm determinant [13] Det (I where {x i , w i } m i=1 is the m-th order quadrature of the interval I such that m i=1 w i f (x i ) computed by the TW system (15) starting from the initial value ϵ = 10 −10 using Mathematica's NDSolve package with WorkingPrecision → 5 MachinePrecision, and that evaluated by the Nyström-type approximation ( 16) with the Gauss-Legendre quadrature of order m = 200, do not exceed 10 −27 for the whole range of variables |a j | ≤ 10, and < 10 −19 for |a j | ≤ 20.Interested readers are invited to verify this statement by running the Notebook Janossy TW N.nb included in Supplementary materials.
Various probability density functions follow from the Jánossy density: The level spacing distribution P (s) [16], the distribution for the nearest neighbor level spacing P nn (t) [20], and the joint distribution for the two consecutive level spacings P c (a 1 , a 2 ) [2] are given by respectively (Fig. 1, right).Finally, the distribution P r (r) for the ratio r = |a 1 |/a 2 of the two consecutive level spacings is given by [2], If we switch the variable from r to r := min(|a 1 |, a 2 )/ max(|a 1 |, a 2 ) = min(r, r −1 ) ∈ [0, 1], its distribution is twice the above P r (r) (Fig. 2), yielding the expectation values for its moments, , 0.4132049292(1), 0.3100223500(1), 0.2460560527(1) 4. Zeroes of the Riemann ζ function.As a showcase of application of our analytic results to the judgement of quantum chaoticity, we compare them to (arguably) the most ideal of all quantum-chaotic spectra: the sequence of zeroes of the Riemann ζ function on the critical line, { 1 2 + iγ n }.Their imaginary parts are supposed to be the eigenvalues of the hypothetical self-adjoint Hilbert-Pólya operator [22].Extensive computational and analytic number theory studies since Montgomery's pair correlation conjecture [23] have fruited in conviction that such an operator, if interpreted as a Hamiltonian, should be ergodic and possess no antiunitary convolutive symmetry, i.e. belong to the unitary universality class of 6/9   [24,25].After unfolding by the Riemann-von Mangoldt formula for the asymptotic density of zeroes, ρ(γ) = 1 2π log γ 2π , the histogram P c (δ − , δ + ) of two unfolded consecutive spacings δ ± = ρ(γ n )(γ n±1 − γ n ) of 10 8 zeroes ending at n = 103800788359 (the largest zero available at the L-functions and modular forms database [26]) perfectly agree with the GUE result (17), as visualized in Fig. 3 (left).Moreover, had we not known the classical formula for ρ(γ) a priori, the perfect match to the GUE could still be deduced from the distributions of the ratios of consecutive spacings of zeroes r n = (γ n+1 − γ n )/(γ n − γ n−1 ) and rn = min(r n , r −1 n ) (Fig. 3, right).Indeed, this observation has been reported in Fig. 4 of the original article [2] that inspired this work.Equipped with the high precision of the moments of r attained by our analytic derivation of P c (a 1 , a 2 ), we can now revisit the rather crude observation of Ref. [2] and improve it to the level of quantifying the systematic convergence of the distribution of the r-ratios of the Riemann zeroes to the GUE result (18).Mean values of the moments of rn in four windows of zeroes [γ N , γ 1.001N +1 ] for N = 10 8 , 10 9 , 10 10 , and 103700788358 are summarized in Table 1.As the height γ N increases, relative deviations from Eq. ( 19), ⟨r k n ⟩/E[r k ] − 1, are indeed observed to vanish systematically from above, heuristically in proportion to ρ(γ N ) −3 , indicating an approach to complete quantum chaoticity in the "thermodynamic limit" γ N → ∞ (Fig. 4).
We remark that the systematic deviations of the statistics of the Riemann zeroes from the GUE (at infinite N ) were previously studied in Refs.[27][28][29].There, the "finite-size corrections" in various statistical distributions for the Riemann zeroes were reported to agree well with those for the circular unitary ensemble (CUE) at finite N eff ≈ 1.446124 • ρ(γ N ).The relationship between these observations and our finding above will need to be clarified.

1 Fig. 4
Fig. 4 Relative deviations of the moments of the r-ratio ⟨r k n ⟩ (k = 1, . . ., 4) in four windows of the Riemann zeroes { 1 2 + iγ n }, n ∈ [N, 1.001N + 1] for N = 10 8 , 10 9 , 10 10 , and ≃ 10 11 , from the GUE results E[r k ], plotted versus ρ(γ N ) −3 .Error bars represent statistical fluctuations due to an arbitrary choice of windows, estimated by the Jackknife method of splitting each data set in 10 bins.Dotted lines are linear fits with the least χ 2 for each set of four points in the same color.

Table 1
Moments of rn for the Riemann zeroes.