Ordered level spacing probability densities

Spectral statistics of quantum systems have been studied in detail using the nearest neighbour level spacings, which for generic chaotic systems follows random matrix theory predictions. In this work, the probability density of the closest neighbour and farther neighbour spacings from a given level are introduced. Analytical predictions are derived using a $3 \times 3$ matrix model. The closest neighbour density is generalized to the $k-$th closest neighbour spacing density, which allows for investigating long-range correlations. For larger $k$ the probability density of $k-$th closest neighbour spacings is well described by a Gaussian. Using these $k-$th closest neighbour spacings we propose the ratio of the closest neighbour to the second closest neighbour as an alternative to the ratio of successive spacings. For a Poissonian spectrum the density of the ratio is flat, whereas for the three Gaussian ensembles repulsion at small values is found. The ordered spacing statistics and their ratio are numerically studied for the integrable circle billiard, the chaotic cardioid billiard, the standard map and the zeroes of the Riemann zeta function. Very good agreement with the predictions is found.


Introduction
Random matrix theory (RMT) describes quite successfully the statistical properties of complex systems spectra of various origins, such as nuclear and condensed matter systems, and has applications in mathematics as well [1,2]. Consecutive level spacing statistics ‡ belong to one of the simplest classes of measures and has been used to differentiate between the spectra of integrable and chaotic origin. Generic integrable systems are expected to follow Poissonian statistics [3], whereas the stastistics of classically strongly chaotic systems are described by the corresponding RMT densities [4]. Other important commonly used fluctuation measures are given by the number variance and spectral rigidity [5][6][7][8]. Their universal behaviours create the opportunity to search for deviations that reveal interesting physics arising in physical systems, e.g. localization or other effects.
In this paper a general class of spacing statistics is introduced, which is based on the ordering of the levels according to the increasing distance from a given level as opposed to consecutive levels. The simplest are the densities of the spacings of the closest neighbour and that of the farther neighbour. One motivation for this comes from applications in the context of perturbation theory [9,10], where it is of importance to know the statistics of the distance of a given level to the closer neighbour (CN) of its two nearest neigbors. The reason is that this gives the strongest contribution due to a smaller energy denominator. The closest neighbour spacing density has already been discussed for the Gaussian unitary ensemble in Refs. [11,12]. The complement of the closest neighbour is the farther neighbour (FN). In general, there is the set of k-th closest neighbour (kCN) and these determine the most significant contributions in the context of perturbation theory and allow for testing longer range correlations in the spectra.
For the closest neighbour and farther neighbour spacing densities random matrix predictions are derived for the Gaussian orthogonal (GOE), Gaussian unitary (GUE), and Gaussian symplectic (GSE) ensembles. The results for Poissonian spectra describing the case of systems with integrable dynamics are also given. These predictions are successfully compared with results for the densities for several disparate systems: the integrable circle billiard, completely chaotic case of the cardioid billiard, the kicked (Floquet) system of the standard map and the statistics of the zeros of the Riemann zeta function, which is known to have properties of the GUE spectra.
One drawback of level spacing densities, especially for many-body systems, is the requirement of unfolding spectra to unit mean spacing. Unfolding works best with an analytical expression for the smoothed density of eigenvalues, such as the Weyl formula for billiards [13], which is not always available. To overcome this difficulty, the ratio of the closer neighbour spacing to the farther neighbour spacing (r = CN/FN) was ‡ Commonly this is referred to as nearest neighbour spacing distribution. However, for the purposes of this paper this would lead to confusion; to avoid this will use consecutive level spacing probability densities throughout the text.
introduced in Ref. [14]. While this ratio has quickly become popular in the context of many body localization as indicator of chaotic behaviour [15][16][17][18][19] and has been studied for various RMT ensembles, the spacings CN and FN themselves are neglected. Based on the closest neighbour and second closest neighbour spacings we propose an alternative, i.e. the statistics of their ratio r CN , defined in Eq. (40), which also does not require an unfolding. As the second closest neighbour could be the farther neighbour or the next-neighbour level in the direction of the closest, this ratio will be in general larger than r.

Ordered spacing probability densities
Consider a sequence of levels {x i } which are assumed to be ordered, i.e. x i+1 ≥ x i , and unfolded [20,4,21]; see Fig. 1. The consecutive level spacings s i = x i+1 − x i have unit average. The following types of ordered spacings will be studied: by determining the k-th element after sorting the set. For k = 1 this reduces to the closest neighbour (CN) and for k = 2 to the second closest neighbour (2CN) spacings. The ordering operation is equivalent to cutting a spectrum at the i th level and reflecting its left or right portion to the opposite side.
Random matrix approximations for the probability density of the farther neighbour and the closest neighbour are derived in sections 2.1 and 2.2 for the GOE, GUE, and GSE, (with Dyson indices β = 1, 2, 4 respectively) based on a 3 × 3 model, and for the Poissonian spectrum. The k-th closest neighbour densities are discussed in section 4. The results for the probability density of the closest neighbour are in agreement with earlier results obtained for the GOE in [12]. In section 3 these results are compared with numerical computations for various dynamical systems and the zeros of the Riemann zeta function. Figure 1. Sketch of the neighbours of x i . In this example x i+1 is the closest neighbour (CN), while x i−1 is the farther neighbour (FN). The second closest neighbour (2CN) is x i+2 , as x i is closer to x i+2 than to x i−1 . Consecutive level spacing probability density is based on the s i = x i+1 − x i .

Farther neighbour (FN) spacings
Start with the joint probability density (JPD) of the consecutive level spacings for the matrix ensemble under consideration. Despite being a highly correlated sequence, the JPD can be written analytically by performing a change of variable in the JPD of the eigenvalues. In principle, all the spacings can be integrated excluding the successive two of interest. An alternative approximate route, in the spirit of Wigner [22], is to consider minimal 3 × 3 random matrix ensembles. Denote the JPD of the three ordered eigenvalues λ 1 ≤ λ 2 ≤ λ 3 by P β (λ 1 , λ 2 , λ 3 ). The FN spacing density is the probability density of max(s 1 , s 2 ) where s 1 = x 2 − x 1 = µ β (λ 2 − λ 1 ) and s 2 = x 3 − x 2 = µ β (λ 3 − λ 2 ) are the two successive spacings scaled by a factor µ β in order to unfold to unit spacing. Such a multiplicative scaling is sufficient as we are dealing only with three levels. The following sequence gives the steps for calculating the FN spacing density for the three Gaussian ensembles: (i) The JPD of three ordered levels {λ 1 , λ 2 , λ 3 } is given by [23] with normalisation N β , and Θ(x) is the Heaviside step function.
(ii) The JPD of the spacings s 1 , s 2 can be calculated from the JPD of the ordered eigenvalues as whereÑ β is a normalisation constant.
(iii) The constantsÑ β and µ β are determined by normalization of P β (s 1 , s 2 ) and requiring the unit average of its marginals, that is Here, P GOE is P 1 etc., to be used where is the explicit abbreviation of the ensemble is favored. Also the scaling factors µ i are related to the constants a, b, c as a = 1/2µ 2 1 , b = 1/µ 2 2 , and c = 2/µ 2 4 . Integrating one spacing results in the consecutive level spacing density p β (s). For example, for the GOE this yields Although different from the Wigner surmise derived from two levels (see Eq. (17) below), the difference is negligible for nearly all applications.
(iv) The cumulative distribution of max (s 1 , s 2 ) is defined as from which the probability density is obtained via where in the final step the symmetry in s i is used. Therefore, the probability densities of the farther neighbour spacings are √ bs e 4bs 2 3 bs 2 − 9 + 3 7bs 2 + 9 , (5b) (vi) For large spacings the behaviour is Thus the small s behaviour is s 3β+1 , while for large s the density is ∼ . In contrast to the high degrees of level repulsion in the Gaussian ensembles, for a Poissonian spectrum it follows using the standard method for finding the probability of the maximum of two random variables, that This goes as ∼ 2s for small spacings, and ∼ 2 exp(−s) for large s, reflecting the larger probabilities of finding both small and large values of FN in comparison to the Gaussian ensembles. The FN behaviour may also be contrasted with the usual behaviour of the consecutive level spacings of the Gaussian ensembles, that is ∼ s β for small s and s β exp(−c β s 2 ) for large s, again indicating the smaller relative variability of the FN.

Closest neighbour (CN) spacings
Following a similar approach, the closest neighbour spacing probability densities can be worked out by first finding the cumulative distribution function for s CN i as Q[min(s 1 , s 2 ) < 0] or alternatively using which implies where again the symmetry of P (s 1 , s 2 ) in Eq. (2) is invoked. From Eq. (4) for the FN density, it follows that where p(s) is the consecutive level spacing probability density. Thus it is an equal mixture of those of the CN and FN.
The probability density functions of CN for the three Gaussian ensembles are For small spacings, the probability densities show level-repulsion, similar to the consecutive level spacing probability densities which goes as ∼ s β , i.e.
A similar calculation yields P CN (s) for a Poissonian spectrum. In general, from the inclusion-exclusion principle it follows that the cumulative distribution function for Z ≡ min(X, Y ) is related to the cumulative distribution functions of X, Y and the joint cumulative distribution of X, Y as, For the Poisson spectrum consistent with the fact that an equal mixture with the FN density in Eq. (8) yields e −s , the first consecutive level spacing density of the Poisson distribution. The result for the Poissonian CN has been used in the context of perturbation theory for the entanglement in bi-partite systems [9,10]. Figure 2 shows a comparison of the analytical results for the closest neighbour densities of Eq. (12) and farther neighbour densities of Eq. (5) with calculations using matrices of the Circular orthogonal ensemble (COE), Circular unitary ensemble (CUE) and Circular symplectic ensemble (CSE); they are generated as described in Ref. [24]. Spectrum Despite approximating the full JPD by a 3 × 3 matrix, the analytical results describe the numerical results very well. This is very similar to how the Wigner surmise captures the consecutive level spacing probability density of the infinite dimensional ensembles from a 2 × 2 matrix ensemble.
The mean values of the closest neighbour and the farther neighbour spacings for the Poissonian, GOE, GUE, and GSE ensembles are summarized in Table 1. Notice that mean values of the two cases are approaching each other with increasing amount of level repulsion in the spectrum, and their sum is always 2, as required by Eq. (11).

Application to various systems
The closest neighbour and farther neighbour densities are examined here for a series of dynamical systems and the Riemann zeta function. For integrable systems the probability density P (s) of consecutive level spacings is given by [3] P Poisson (s) = e −s .
In this case the density does not vanish at the origin, i.e. P (s) → 1 for s → 0. In contrast, the densities of classically strongly chaotic systems follow the corresponding RMT densities [4]. The GOE result for the consecutive level spacing density is to a very good approximation given by For the GUE one has to a very good approximation For these two cases and the GSE the spectrum displays level repulsion as P (s) ∼ s β for small s.

Quantum billiards
A classical billiard system is given by the free motion of a point particle inside some domain Ω with specular reflections at the boundary, i.e. the angle of incidence equals the angle of reflection. Depending on the shape of the boundary, the dynamics can range from integrable such as for the circular, rectangular or elliptical billiards, to fully chaotic motion as for the Sinai-biliard [25], Bunimovich stadium billiard [26], or the cardioid billiard [27][28][29][30].
A corresponding quantum billiard is described by the stationary Schrödinger equation (in units = 2m = 1) − ∆ψ n (q) = E n ψ n (q) , q ∈ Ω, with Dirichlet boundary conditions, i.e. ψ n (q) = 0 for q ∈ ∂Ω. The unfolded spectrum is obtained by the mapping where N (E) is the mean behaviour of the spectral staircase function and is given by the generalized Weyl formula [13] N where A denotes the area of the billiard, and L := L − − L + , where L − and L + are the lengths of the pieces of the boundary ∂Ω with Dirichlet and Neumann boundary conditions, respectively. The constant C takes curvature and corner corrections into account.
As an example integrable system consider the circular billiard with unit radius and Dirichlet boundary conditions. For the desymmetrized (half circular) billiard and Neumann boundary conditions on the symmetry line, the eigenvalues are E kl = j 2 kl , where j kl is the l-th zero of the Bessel function J k (x) with k = 0, 1, 2, ... and l = 1, 2, .... The corresponding constants in Eq. (21) are given by A = π 2 , L = π − 2 and C = − 1 24 . The first 10 5 eigenvalues give for the closest neighbour spacing density P CN (s) and  farther neighbour spacing density P FN (s) the results shown in Fig. 3. Good agreement is found with the predictions of Eqs. (15) and (8), respectively. Next consider the chaotic cardioid billiard, which is the limiting case of a billiard family introduced in Ref. [27]. Its boundary is given in polar coordinates (ρ, φ) by The cardioid billiard is proven to be strongly chaotic [28,[31][32][33][34]. For the desymmetrized cardioid billiard with Dirichlet boundary conditions, the statistics is based on the first 11,000 eigenvalues [35], determined using the conformal mapping technique [36,37]. For a detailed investigation of the spectral statistics see [38]. For the mean behaviour N (E), Eq. (21), of the spectral staircase function the constants are A = 3π/4, L = 6, and C = 3/16. Figure 4 shows the closest neighbour spacing density P CN (s) and farther neighbour spacing density P FN (s) and there is good agreement with the predictions of Eqs. (12a) and (5a), respectively. The comparison with the consecutive level spacing density P GOE (s) shows that P CN (s) has its maximum for smaller s and P FN (s) has its maximum for larger s and the mean values, as shown in Table 1, are also below and above 1.

Quantum maps
Another extensively studied class of systems are quantum maps which are unitary operators on a finite-dimensional Hilbert space. Their classical limit is given by a discrete-time mapping on phase space. A prototypical example is the so-called kickedrotor [39]. It is obtained from the Hamiltonian where the sum describes a periodic sequence of kicks with unit time as the kicking period and V (q) = K cos(2πq)/4π 2 . Considering the dynamics immediately prior to consecutive kicks leads to the standard map, which is a two-dimensional area-preserving map. The map is considered on a twodimensional torus with periodic boundary conditions in both q and p. If the kicking strength K is sufficiently large, the map is strongly chaotic with a Lyapunov exponent of ≈ ln(K/2) [39]. This excludes some special values corresponding to small scale islands formed by the so-called accelerator modes [39]. This is supported by recent rigorous results showing that the stochastic sea of the standard map has full Hausdorff dimension for sufficiently large generic parameters [40]. As example we use K = 10 for which the dynamics of the map numerically looks chaotic, i.e. no regular islands are observed on any relevant scales. The quantum mechanics on a torus phase space leads to a finite Hilbert space of dimension N , see e.g. Refs. [41][42][43][44][45]. The effective Planck constant is h = 1/N . The quantum map U is an unitary operator, where n, n ∈ {0, 1, ..., N − 1}. The parameters β and α are the quantum phases due to periodicity in position and momentum respectively. The cases when these are 0 or 1/2 correspond to periodic and anti-periodic boundary conditions respectively. We choose (α, β) = (0.2, 0.3) ensuring that parity and time reversal symmetry are both broken.
Thus, the spectral statistics should follow those of the circular unitary ensemble (same as GUE). Figure 5 shows the closest neighbour spacing density P CN (s) and the farther neighbour spacing density P FN (s) for the quantized standard map for N = 30, 000. The agreement with the expected GUE results, Eqs. (12b) and (5b), is very good.

Riemann zeta function
An interesting connection between the spectral properties of chaotic dynamical systems and number theory arises with the statistics of the zeros of the Riemann zeta function. According to a famous conjecture, its zeros are expected to lie on the critical line [46], i.e. ζ(1/2 + it n ) = 0 .
Their counting function is given by the Riemann-von Mangoldt function [46][47][48] As is well known for the spacing densities, the convergence to the CUE results is rather slow [49]. Therefore, the 10, 000 zeros starting from the (10 22 + 1)-th [50] are used. To numerically perform the unfolding for the tabulated zeros, which have an offset δ = 1370919909931995300000, this has to be accounted for when applying (27).
where the approximation ln(1 + x/δ) ≈ x/δ has been used. After this unfolding the mean spacing is 1.00011. Figure 6 shows the closest neighbour P CN (s) and the farther neighbour spacing densities P FN (s) and very good agreement with the expected GUE results, Eq. (12b) and Eq. (5b), respectively, is found.

k-th closest neighbour (kCN) spacings
By means of the k-th closest neighbour spacing density longer-range statistical properties can be investigated.

Poisson sequences
For the case of a Poissonian spectrum it is possible to obtain an analytical prediction. The steps of the proof are illustrated for the 2CN: • The unfolded eigenvalues form a unit mean Poisson process such that the consecutive level spacings s i are distributed according to exp(−s). • Each semi-infinite sequence obtained by splitting this infinite sequence of levels from a fixed level is also described by a unit mean Poisson process. • Flipping one of these two semi-infinite sequences and merging results in a semiinfinite sequence described by a Poisson process with rescaled spacing, 1/2. • Therefore, the consecutive level spacings for this new semi-infinite sequence would be the same as the original infinite sequence, except for the change of the mean spacing. Hence, for example exp(−s) → 2 exp(−2s). Explicitly, The means are given by µ Poisson k = s kCN = k/2 (30) and the variances by Figure 7 shows P kCN (s) for k = 2, 3, 4, 5 for the integrable circle billiard in comparison with Eq. (29). Overall, very good agreement is found, but with increasing k some very small deviations begin to show. They are a manifestation of the increasing deviations of the variance for larger k from k/4; see Fig. 8. This is similar to the well-known saturation of the spectral rigidity and number variance [7,8,51], which also determines the variance of the k-th consecutive level spacing statistics [6].  Figure 9 shows the k-th closest neighbour spacing density for the cardioid billiard. From its definition, it is clear that the likelihood of the distance to the k-th neighbour is the sum of correlated, yet randomly behaving consecutive level spacings. Supposing a central limit theorem for increasing k, the density is expected to become more and more Gaussian [6] G(k, s)

RMT sequences
This is in very good agreement with the numerically obtained histograms. In addition, it is possible to give heuristic arguments for the means and variances of the kCN. Even though correlations exist in the spectra of chaotic systems, the definition of kCN is equivalent to cutting an infinite sequence at a level and flipping one of the two semi-infinite subsequences just as for the Poisson case. The mean spacing reduces to 1/2 as before. However, the spectral correlations introduce some changes. First, relative to the result for the Poisson sequence of Eq. (30), there is a shift of the mean spacing. RMT spectra exhibit a correlation hole. Consider the case before superposing the flipped sequence. Due to the symmetric, nearly evenly spaced sequence of P kCN (s), the mean of P kCN (s) is approximately equal to the width s of the interval required to find k spacings on average. This gives implicitly After superposing the left and right sequences, the density is doubled, and the mean is half of the result. Rearranging terms, this gives An approximate form of which follows by replacing the upper integration limit with the leading dependence of s (≈ k/2), Using exact expressions of Y 2 (s) [6, Eq. (L1)] for the respective ensembles, we obtain These expressions are compared with results from the cardioid billiard and standard map, respectively, in Fig. 10(a). One sees that the trend towards k/2 − 1/4 is confirmed.
The behaviour for small k is fairly good, but shows some slight deviations. In particular for the standard map, an odd-even effect is more pronounced than for the theoretical prediction.
The scale of the variances can be given by noting their similarity to the number variance. Neglecting three-point correlations and the sawtooth corrections [52], the spacing variances should be roughly a quarter of the number variance of appropriate argument. Thus, they should grow no more than logarithmically with increasing k. For an illustration of these quantities and their comparison to these approximate forms, see Fig. 10(b). Considering the simplifying assumptions just invoked, for the quantized standard map the variance σ 2 k follows 1 4 Σ 2 GUE ((k + 1)/2) up to a rather small shift very well. However, the agreement of σ 2 k for the cardioid billiard with 1 4 Σ 2 GOE ((k + 1)/2) is not quite as good. The deviations for larger k correspond to the early saturation of the number variance observed for this system [38,53]. To emphasize the distinction, the figure also shows results obtained for an ensemble of 500 COE matrices of dimension 100, which again agrees very well with the prediction up to a small shift. To give an idea of the role of 3-point correlations in the P kCN (s), the example of the second closest neighbour spacing density is calculated with the approximation that the spacings themselves are independent. This incorporates most of the 2-point correlation information, but completely removes 3-point correlations. By definition no more than two levels from either side of a chosen level can be the second closest neighbour, thus isolating 5 levels {x i−2 , x i−1 , x i , x i+1 , x i+2 } of the spectrum; see Fig. 1.
spacings has two elements as the consecutive level spacings, namely x i − x i−1 = s i−1 and x i+1 − x i = s i , and two as sum of the consecutive level spacings, namely The second closest neighbour spacing density can be obtained by calculating the order statistics of the correlated spacings defined earlier. To a first approximation, this sequence can be treated as being formed by independent, but not identical random variables. Consecutive level spacings themselves are distributed according to the Wigner surmise, and sum of two are calculated by convolution. These approximations give a form for the density and the distribution (cumulative density) of the elements of the sequence formed by consecutive level spacings and their sum. Using the Bapat-Beg theorem [54] for the density of k th order statistics in terms of the permanent of a matrix constructed using the probability densities and cumulative distribution functions of the individual elements, the 2CN spacing density is, These approximate densities in the limit of small and large s behaves as, The approximation, Eq. (37), is also shown in Fig. 9(a) and some significant deviation is found. Nevertheless, the approximation reproduces the mean 2CN within ∼ 3%. Similar calculations can be carried out for the GUE ensemble, leading to a rather lengthy expression, which is omitted here. The quality of agreement with the numerical results is similar to the GOE case.

Ratio statistics of CN/2CN
To avoid the unfolding procedure, the ratio statistics of successive level spacings has been introduced as [14] which in terms of the ordered spacings is the same as the ratio s CN i /s FN i . Notice that this ratio mostly contains information of up to 3-point correlations. It has become very widely used in the context of many-body systems, where the unfolding may be difficult to carry out properly.
An alternative to the above measure with a cleaner interpretation in terms of the ratio of the spacing of the two closest levels is given by As for the ratio statistics, Eq. (39), this quantity does not require an unfolding of the levels. There are circumstances in a spectrum for which the farther neighbour would be the second closest neighbour (especially for the GUE and GSE), and r CN i coincides with r i , but this does not always occur. The distinctions between these closely related measures are nicely reflected in the different densities, FN and 2CN, of the previous section. Though an analytical prediction for this quantity is not given, some salient features can be quickly inferred from this interpretation. For example, the occurrence of eigenvalues in integrable systems is described by a Poisson process which means that the probability of finding a level is independent of the presence or absence of another level, therefore the probability of the ratio of r CN i would be equal for all ratios. This is clearly seen in Fig. 11 (a). For the RMT ensembles, level repulsion implies zero probability for the ratios r CN at the origin. To lowest order, the behaviour near small   Table 2. Numerical mean values r and r CN of CN/2CN and CN/FN, respectively. ratios is P (r CN ) ∝ (r CN ) β . These features are visible in the numerical histograms of Fig. 11 (b), (c) and (d). The comparison of the mean ratios of r and r CN for the circle billiard, cardioid billiard, standard map, and Riemann zeroes are given in Table 2. The increase of the mean value of P (r CN ) relative to P (r) is consistent with the fact that the mean of P 2CN (s) is smaller than P FN (s) and both of them are continuous, smooth densities.

Summary and discussion
In this paper, the closest neighbour and farther neighbour spacing probability densities are introduced and predictions obtained for the random matrix ensembles GOE, GUE, and GSE based on a 3 × 3 matrix modeling, and Poisson spectra. The results are successfully compared with numerical computations for the integrable circle billiard, the fully chaotic cardioid billiard, the standard map with chaotic dynamics and broken time reversal symmetry, and the zeros of the Riemann zeta function. More generally, the class of k-th closest neighbour spacing statistics are discussed, which are numerically found to be well described by Gaussians with mean increasing by half and variances depending on the case. One application of the closest neighbour spacing density is in the context of perturbation theory, and the farther neighbour sheds light on understanding the ratio defined using successive level spacings.
In addition, the statistics of the ratio of the closest neighbour spacing to that of the second closest neighbour is introduced as alternative measure for the correlation properties of complex quantum spectra, which does not require unfolding. This ratio has a simpler interpretation. Finding an analytical prediction for the probability density of the ratio is an interesting task left for the future. A straight-forward generalization, left for future work, would be to consider the statistics of the ratios s kCN i /s (k+1)CN i , which would give insight into longer range correlations. It is noteworthy that a recent work [55] explores certain higher range spacing ratios of a kind that is not based on ordering.