Density and spacings for the energy levels of quadratic Fermi operators

The work presents a proof of convergence of the density of energy levels to a Gaussian distribution for a wide class of quadratic forms of Fermi operators. This general result applies also to quadratic operators with disorder, e.g., containing random coefficients. The spacing distribution of the unfolded spectrum is investigated numerically. For generic systems the level spacings behave as the spacings in a Poisson process. Level clustering persists in presence of disorder.


I. INTRODUCTION
In a variety of situations, one encounters quadratic forms in Fermi operators where the Fermi operators c i 's obey the canonical anticommutation relations {c i , c j } = 0 and {c i , c j † } = δ ij , and the coefficients satisfy A ij = A ji ∈ R and B ij = −B ij ∈ R for i, j = 1, 2, . . .. The quadratic form (1) defines a symmetric operator H n acting on a Hilbert space of dimension 2 n . This operator represents the Hamiltonian of a system of quasifree fermions.
Such quadratic operators are of fundamental interest for several reasons. First and foremost, these operators can be diagonalized exactly using an explicit normal mode decomposition, 23 and certain quadratic Hamiltonians are good approximations for more complicated two-body interactions. Furthermore, fermionic models share a close relationship to interacting spins in one dimension, and they are among the simplest systems in which quantum phase transitions occur and entanglement measures can be computed. The literature on quasifree fermions, their relationship to spin systems, and to other areas of physics, such as conformal field theory and random matrix theory, is immensely vast. For a review, see Refs. 1, 10, 12, and 20. In the last decades, models with disorder, i.e., containing random parameters, have also been considered. 15,28,30,31 For quadratic forms in Fermi operators, it is of interest to know whether, in the limit of large n, the density of energy levels and the spacing distribution converge to identify the limit. Of course, the spacing distribution of the unfolded spectrum requires the knowledge of the density of energy levels. We discuss these questions for a broad class of systems that includes the following examples: This representation is exact in the case of a spin chain with free ends. (ii) A quantum bond percolation model on a lattice Γ with n sites consisting of a tight-bind Hamiltonian of the form 7,33,34 where the summation runs over nearest neighbour sites, and the hopping matrix elements t ij ∈ R are independent Bernoulli random variables Pr(t ij = 1) = 1 − Pr(t ij = 0) = p ∈ (0, 1). (iii) The Anderson model 2 is one of the simplest models incorporating the essential competition between the hopping term (discrete Laplacian) and the on-site disorder (random potential). For a generic lattice Γ, the Anderson Hamiltonian for noninteracting fermions can be written as with random on-site potential v i ∈ R; usually v i 's are independent with mean zero and finite variance W 2 . (iv) More general non-sparse random quasifree fermions Hamiltonian. For instance, one can consider the Hamiltonian (1) with A ij and B ij independent Gaussian variables, modulo the symmetries A ij = A ji and B ij = B ji . It turns out that this model is related to the real Ginibre ensemble of random matrix theory. 18 In the traditional paradigm of condensed matter physics, the number of particles is so large that questions on the macroscopic density of energy levels, i.e., the behaviour of the energy levels in the "bulk" very far from the ground state, are meaningless. The situation has changed recently. Over the past few years, experimental developments have allowed the study of systems with a small and controlled number of particles, and therefore, a direct measure of the level density might be within the reach of current experimental capabilities.
These considerations have recently triggered the attention of some authors on the problem of convergence and universality of the limiting level density of many body systems. In two pioneering papers, Hartmann, Mahler, and Hess 19 considered generic many body quantum systems with nearest neighbour interaction. They proved that, provided that the energy per particle has an upper bound, the energy distribution for almost every product state becomes a Gaussian in the limit of infinite number of particles. More recently, Atas and Bogomolny 3,4 investigated numerically and theoretically the energy levels of several interacting spin 1/2 systems and concluded that the density of levels converges to a Gaussian. Using an adaptation of the line of reasoning in Ref. 19,Keating,Linden,and Wells 21,22,37 proved convergence to a Gaussian distribution for spin chains with generic pair interactions, including the case of spin glasses, i.e., interaction with random couplings. This result has been extended to spin systems on more general graphs by Erdös and Schröder. 13 The algebraic identities satisfied by the Pauli matrices representing spin 1/2's play a key role in the proofs in Refs. 13, 21, 22, and 37. Our goal here is to show that the density of energy levels of a wide class of quadratic Fermi operators (both deterministic and random) converges to a Gaussian distribution in the limit of large n. The proof of this universal result relies on the connection between the spectrum of H n and the subset-sum structure arising in the normal mode decomposition. This result explains some of the previous conjectural statements and numerical observations by Atas and Bogomolny on the level density of certain (nonrandom) spin systems. Additionally we provide a uniform bound (based on a Berry-Esseen inequality) on the rate of convergence.
We also consider the level spacing distribution of such operators. Numerical investigation shows that both deterministic and random models exhibit level clustering (Poisson statistics); this behavior is compatible with the celebrated Berry-Tabor philosophy for generic integrable systems, 8 even in presence of disorder. In the course of the paper, we also present a few short examples illustrating the general theorems.
The paper is organised as follows. In Sec. II, we set the notation and review the consequences of the normal mode decomposition. Then, in Sec. III, we present our main results on the limiting density of energy levels and the rate of convergence to the limit. In Secs. IV-VII, we apply the general theorems to the examples (i), (ii), (iii), and (iv) discussed above, thus illustrating in physical models the universality of Theorem 1 and Corollary 1. Finally, in Sec. VIII, we present the numerical observations on the level spacing distribution.

II. GENERALITIES ON FERMI OPERATORS
Let us order the 2 n eigenvalues of H n as where the normal modes η k , η † k are the Fermi operators, and the elementary excitations λ k,n ≥ 0 are the singular values of P n (A + B)P n and K n = TrP n AP n /2.
The following well-known properties of the Fermi operators η k and η † k are immediate consequences of the canonical anticommutation relations. 27 First, the η † k η k are Hermitian operators with eigenvalues 0 and 1. Second, η k (η † k ) acts as a lowering (raising) operator on the normalised eigenvectors of η † k η k with eigenvalue 1 (0). Moreover, the η † k η k 's form a set of mutually commuting operators and therefore they can be simultaneously diagonalized. These three facts imply that there exists a normalised vector |0 (the vacuum state) which is an eigenvector of all the η † k η k 's with corresponding zero eigenvalue: η † k η k |0 = 0. A set of 2 n normalised eigenvectors of η † k η k (k = 1, . . . , n) can be built up by exciting the vacuum state; the normalised vector |α 1 α 2 . . . α n = (η † 1 ) α 1 (η † 2 ) α 2 · · · (η † n ) α n |0 , with α k = 0 or 1, is an eigenvector of η † k η k with eigenvalue α k . Therefore, from (6) we have The spectrum of H n is constructed by exciting the ground state energy E 1,n = K n − 1/2 k λ k,n by the elementary excitations λ k,n . Hence the spectrum is characterised in terms of the subset sums of elementary excitations as follows: E is an eigenvalue of H n if and only if The density of energy levels is defined as the empirical normalised measure and from (8) it follows that Up to a shift, the empirical measure of E k ,n is given by the distribution of the sum of n independent variables r 1 λ 1,n , . . . , r n λ n,n . In fact, it is possible to compute the Fourier transform of (10) This computation shows that the empirical distribution of the energy levels E k is the distribution of a sum of independent random variables. It is then plausible that for large n, after a suitable rescaling, the distribution of energy levels converges to a Gaussian. After all, the many body Hamiltonian (6) is a sum of single particle (commuting) operators, and the total spectrum is given by the sum of the individual spectra. In Sec. III, we specify exact conditions for this convergence. Note that the variables r k λ k 's are independent but not identically distributed, e.g., E[r k λ k,n ] = 0 and E[(r k λ k ) 2 ] = λ 2 k,n /4. Before stating the main theorems we conclude this section with a last computation to prepare the ground to what follows. If we knew that the limiting level density is Gaussian then the limit would be identified by its mean and variance. The moments of the counting measure (9) are related to traces of powers of H n by the following identity: In particular, mean and variance are given by the traces of the first two powers TrH n and TrH 2 n . A direct computation of these traces is possible using Wick's calculus. The only non-traceless products of Fermi operators that we need are Using (13)-(15) one finds The above quantities are mean and variance of the finite-n level density.

Theorem 1 (Density of energy levels).
Let H n be the quadratic form (1). Assume that, denoting X n = P n (A + B)P n , the following conditions are true: Tr(X T n X n ) = σ 2 < ∞.
Then, the density of shifted and rescaled energy levels weakly converges, as n → ∞, to a centred Gaussian probability measure with variance σ 2 , Theorem 1 can be proved by checking the Feller-Lindeberg conditions 14 in the central limit theorem for independent nonidentical random variables. We present however a more direct proof based on elementary computations. Hypothesis (i) of Theorem 1 can be rephrased as λ k,n = o(n 1/4 ) for all k, (19) meaning that the elementary excitations do not grow too fast with n. This assumption is similar (but in a weaker sense) to the condition of finite energy per particle in Hartmann, Mahler, and Hess theorem. 19 Note also that according to (16). Hypothesis (ii) is thus a condition on the second moment of the density of energy levels.
Proof of Theorem 1. Let λ k,n (k = 1, . . . , n) be the singular values of X n . Repeating the computation in (11) we find The key point to appraise (21) is the following identity.
Claim. Let (u i ) i∈N be a sequence of complex numbers such that and the following limit exists and is finite: Then (A generalization of the identity lim If we accept the claim, we can prove the theorem as follows. For any fixed t ∈ R n k=1 cos tλ k, By the claim and hypotheses (i) and (ii) (using max k λ 2 k,n = X T n X n op ), the last expression converges to exp(−σ 2 t 2 /2), and by Lévy's continuity theorem, this proves (18).
It remains to prove the claim. We adapt a proof given in [Ref. 24, Lemma A.5]. Set Note that the function log(1 + z) + z has a double zero at z = 0. Hence is analytic in the open disk |z| < 1 (in particular it is continuous and bounded). The finite product can be written as Therefore By continuity, there exists 0 < R < 1 such that L(R) = 1. From (22) it follows that, for n sufficiently large, M n /n ≤ R, and by the maximum principle, |L( u i n )| ≤ 1. We conclude that, for large n since for any z, |e z − 1| ≤ |z|e |z | . By (22) and (23) the above inequality implies the claim (24).
The following result provides a uniform bound on the discrepancy between the finite-N empirical density of energy levels and the limiting Gaussian (in the sense of the Kolmogorov distance between probability distributions).

Proposition 1. Denote
Then for all n the following bound on the distance between the counting measure of the normalised energy levels E k ,n /s n and the standard Gaussian distribution holds for an absolute constant C that may be chosen as C = 6.
Proof. To prove (32), we use a classical Berry-Esseen inequality for independent nonidentically distributed variables. The empirical distribution of the shifted energy levels is the same as the empirical distribution of the sum of independent centred random variables x 1 , . . . , x n with the position where r 1 , r 2 , . . . are i.i.d. binary variables (note that the x k 's are not identically distributed). Denote by F n the cumulative distribution of the normalised sum (x 1 + · · · + x n )/( n k=1 E[x 2 k ]) 1/2 . Then for all x and n where C ≤ 6 (see Ch. XVI.5, Theorem 2 in Ref. 14). An elementary computation shows that where the expectation value is taken with respect to r 1 , . . . , r n . This concludes the proof since the λ k,n 's are (up to a rescaling) the singular values of X n .
Example 1. We show that the rate n 1/2 in (32) is optimal. Suppose that A ij = ξδ ij (ξ ∈ R) and B ij = 0. Hence, the quadratic form reads In this case, the elementary excitations λ k,n are all equal to ξ and the empirical distribution of the energy level E k ,n is given by the distribution of the sum of i.i.d. variables x k = ξr k with r k as above. Therefore we have  (32).
Theorem 1 can be adapted to deal with random quadratic Fermi Hamiltonians (see the examples (ii), (iii), and (iv) presented in the Introduction). Let (Ω, F, P) be a probability space. The expectation with respect to P will be denoted by E. Let us suppose that A(ω) and B(ω) (ω ∈ Ω) are random double arrays of real numbers satisfying A(ω) ij = A(ω) ji and B(ω) ij = −B(ω) ji . Hence (1) defines a sequence of random quadratic forms H n (ω) in Fermi operators. Our approach to proving convergence to a Gaussian limit consists of two steps: first, we average over fictitious binary variables (using Theorem 1) for a given realization of the disorder (A(ω) ij and B(ω) ij ); then, if the first average in the limit of large n is independent of the realization ω, we can average over the disorder (i.e., with respect to P). Note that all random variables are defined on the same probability space. We have the following result as a corollary of Theorem 1.
Then, the sequence of the density of rescaled energy levels weakly converges in average, as n → ∞, to a centred Gaussian probability measure with variance σ 2 . (This means that as n → ∞, for all f bounded and continuous.) Proof. The proof is based on the representation of the shifted energy levels in terms of the set of fictitious independent binary variables r k E k,n (ω) − K n (ω) = n k=1 r k λ k,n (ω), where λ k,n (ω) are the singular values of X n (ω). Let us introduce the sets By hypothesis, P(S i ) = 1 for i = 1, 2, and therefore P(S) = 1. Hence, if ω ∈ S, by Theorem 1 The above convergence holds P-almost surely (for all ω ∈ S). Moreover, the function x → exp(ix) is absolutely bounded and therefore the almost sure convergence can be promoted to convergence in mean The proof is completed by using Lévy's continuity theorem.
Classes of random matrix ensembles which include quantum spin glasses (random two-spin interaction) on generic graphs have been recently considered in Refs. 13, 21, 22, and 37. For these Hamiltonians, using the algebraic identities for Pauli matrices, it has been proved that the limiting spectral density, as the graph cardinality increases, is Gaussian. For spin 1/2's with nearest neighbourhood random interaction, by the Jordan-Wigner transformation, those systems are equivalent to random quadratic forms of Fermi operators and our method provides an alternative proof of these results. For generic Hamiltonians H n , the inverse Jordan-Wigner transformation maps the problem to spin 1/2 systems with more complicated interactions not considered in previous works. Our method of proof relies on the subset sum structure in the spectrum of quadratic Fermi operators, and it is sufficiently robust to be extended to a large class of random Hamiltonians. An exceptional example of a random Hamiltonian that does not exhibit a Gaussian limit is presented below.
In Secs. IV-VII, we discuss explicit examples in detail. In Sec. IV, we illustrate the method on spin 1/2's systems with fixed nonrandom couplings. We consider in detail the XY model with free boundary conditions and the Ising model with the transverse field studied in Ref. 3. Then, we present our results for the quantum percolation models and the Anderson models (Sec. V). In Sec. VI, we establish the connection between non-sparse Gaussian quadratic operators and the Ginibre ensemble of random matrix theory. Finally, in Sec. VII, we apply our theorem to other random band models. The level spacing distribution is discussed in Sec. VIII.

IV. SPIN CHAINS
As described in the Introduction, chains of interacting spin 1/2's can be mapped to systems of spinless fermions. We shall apply our theorems to those systems.
The paradigmatic example is provided by the XY chain, a canonical toy model for quantum spin systems routinely used as a first example to illustrate new concepts. Assuming free boundary conditions, the Hamiltonian of the XY -model for n spins can be written as (2). In this model, A and B have a tridiagonal form Note that A ii = 0 (hence K n = 0). The elementary excitations are 23 and therefore, the density of energy levels of the XY model converges to Similar considerations can be extended in the presence of external fields. For simplicity we consider the Ising model in the transverse field − n j=1 σ x j σ x j+1 + hσ z j , where h ≥ 0 is the external magnetic field. The problem can be reduced to a quadratic form in Fermi operators whose normal mode decomposition has with phases θ k,n = 2π(k−1) n − π, (k = 1, . . . , n) equidistributed. Now |λ k,n | ≤ (2 + h) and We conclude that according to [Ref. 3,Eq. (20)]. Of course, at zero magnetic field h = 0, we recover the limit density (48) of the XY model in the Ising limit γ → 1.

A quantum bond percolation model (3) can be cast in the form
where i, j ∈ V denotes the vertices (sites) of a graph Γ = (V , E), and A ij (ω) = t ij (ω)1 (i,j)∈E is the adjacency matrix of Γ weighted by random independent Bernoulli variables P(t ij = 1) = 1 − P(t ij = 0) = p ∈ (0, 1) on a probability space (Ω, F, P). We assume that the graph Γ is a connected regular lattice; in particular, Γ does not contain loops (therefore A ii = 0) and the degree of the vertices is constant (d is called the coordination number of the lattice.) It is well known that the largest eigenvalue of the adjacency matrix of a graph is bounded by the maximal degree. This implies that X T n X n op ≤ d. We then compute for P-almost all ω. Therefore, by Corollary 1, we have for P-almost all ω.

VI. GAUSSIAN QUADRATIC FORMS AND THE GINIBRE ENSEMBLE
Let us consider the Hamiltonian (1) with random coefficients A ij (ω), B ij (ω), ω ∈ Ω. We consider the case of A ij (ω) = a ij (ω)/ √ n and B ij (ω) = b ij (ω)/ √ n independent Gaussian variables, modulo the symmetries a ij = a ji and b ij = b ji with mean and variance For Gaussian random variables, the problem is simplified thanks to the following observation: if Z 1 , Z 2 are independent and identically distributed normal variables, then (Z 1 + Z 2 ) and (Z 1 Z 2 ) are independent normal variables. Therefore the entries of the n × n matrix X n (ω) = P n (A(ω) + B(ω))P n are i.i.d. Gaussian variables; hence where G ij are the i.i.d. standard real Gaussian variable. (G is a random matrix belonging to the real Ginibre ensemble. 18 ) We have therefore established that the elementary excitations λ k,n (k = 1, . . . , n) of a quadratic form with i.i.d. Gaussian coefficients are distributed as the singular values of a real Ginibre matrix G(ω) of size n (equivalently, λ 2 k,n are the eigenvalues of a real n × n Wishart matrix W(ω) = G T (ω)G(ω)).
It is well-known that the singular values of n × n Ginibre matrix whose entries are O(1) are typically of order O( √ n). We therefore rescaled the coefficients A ij , B ij by √ n to get a sensible limit for the density of energy levels. In fact, using classical asymptotic results on the extreme singular values of random matrices with i.i.d. entries, 5,17 we know that with probability 1 all the elementary excitations λ k,n (ω) lie in a fixed interval for large n. More precisely we have for P-almost all ω. By the strong law of large numbers, we also have lim n→∞ 1 4n Tr(X T n (ω)X n (ω)) = lim n→∞ s 2 2n 2 n i,j=1 for P-almost all ω, and by Corollary 1, we conclude that for n → ∞ In the rest of this section, we use the relation with the Ginibre ensemble to obtain results on the ground state energy and energy gap. The steps of proof are elementary and they borrow the difficult technical statements from previously known results in random matrix theory. 26 Under the above assumptions of A ij and B ij , we have that n 2s 2 (λ 1,n , λ 2,n , . . . , λ n,n ) where the joint probability density of the n random variables x k 's is The joint law (62) From these results we derive now a few properties of the ground state of H n . The ground state energy is the lowest level E 1,n and we denote by ∆ n = E 2,n − E 1,n the ground state energy gap.
(ii) The rescaled energy gap n/2s 2 ∆ n converges in distribution to a random variable whose probability density function is By the same proof one shows the almost sure convergence of the rescaled largest energy level n 3/2 E 2 n ,n . Therefore, the numerical range of H n is roughly (an 3/2 , an 3/2 ) with a = (2 3/2 /3π)s. Note that ∆ n = O(n −1/2 ). Hence the system is gapless.
Proof of Proposition 2. The ground state energy is given by (see (8)) By the law of large numbers n −3/2 K n = n −2 n i=1 a ii converges to zero almost surely. Using (61) and the quarter law (63) the following almost sure convergence holds: This proves (64). The ground state energy gap is given by the smallest elementary excitation Remark. We expect that one could generalize this analysis to non-Gaussian variables whose first four moments match the Gaussian moments using Lindenberg exchange strategy. Using this technique one can replace the Gaussian variables a ij and b ij one at a time by random variables from a desired distribution. This approach is widely used to prove versions of the four moment theorem. 35

VII. OTHER RANDOM BAND QUADRATIC FORMS
In this section, we show that Corollary 1 applies to the case when A and B are random band arrays. We introduce a parameter W n ≥ 1 which corresponds to the number of non-zero diagonals, i.e., A ij = B ij = 0 if |i − j| > W n .
We normalize A ij (ω) and B ij (ω) to ensure that condition (ii) of Corollary 1 is satisfied. To compute the normalization of the matrix entries in terms of W n , we want Tr(X T n (ω)X n (ω)) = lim n→∞ 1 4n to be finite and non-random. If there are W n non-zero diagonals, the matrix X n has on the order of nW n non-zero entries, in the sense that we can take A ij (ω) = a ij (ω)/ √ W n and B ij (ω) = b ij (ω)/ √ W n with a ij and b ij i.i.d. standardized random variables to achieve the finite limit in (69). Here we do not need a ij and b ij to be Gaussian.
We now show that condition (i) of Corollary 1 is also satisfied and therefore the density of energy levels of random band quadratic forms converges to a Gaussian. Suppose that W n = o(n 1/2 ) and that a ij (and b ij ) has exponential decay, in the sense that there exists δ > 0 such that Ee δ |a ij | < ∞. Then, letting X n = P n (A + B)P n , lim n→∞ n −1/4 X n op = 0 P-a.s. .
We proceed to a proof of (70) by showing that for all L > 0, which implies (70) by the Borel-Cantelli lemma. Note that by triangle inequality, X n op ≤ P n AP n op + P n BP n op . The argument will be identical for the two terms on the right hand side so we will focus on the first one. For a symmetric matrix, the operator norm is equal to the largest modulus of the eigenvalues and it is therefore dominated by any matrix norm. In particular P n AP n op = sup Then using that the Z i 's are identically distributed, by the union bound, we obtain Since a 1j , j = 1, . . . , n, are i.i.d. random variables, we can apply a Chernoff bound to get by which we conclude that (71) holds true. For more general sharp concentration inequalities on the operator norm of random matrices see, for instance, Ref. 6.

VIII. LEVEL CLUSTERING
One of the most commonly studied statistical measures of a given spectrum is the level spacing distribution P(x), i.e., the distribution of gaps between consecutive levels. The first step to unravel meaningful information from the spacings is to unfold the spectrum in such a way that the average level spacing in the neighbourhood of each transformed level is unity. In other words, the unfolding procedure is the scaling transformation that removes the irrelevant effects of the varying local mean density. A natural way to unfold the spectrum is by mapping each level E k ,n into a new variable e k ,n defined as the fraction of energy levels in the spectrum below E k ,n . In practice, the variation of the density of levels needed for the unfolding is included by fitting the integrated level density, or, when explicitly known, by using the limiting level density as an approximation. We have numerically studied the level spacing distribution for a few instances of quadratic Fermi operators.  conditions. As illustrated in the left panel, the histogram representing the numerical empirical measure of the energy levels is almost indistinguishable from the limiting Gaussian density. For this reason, we have used the limiting Gaussian density in (48) to unfold the spectrum. P(x) of the unfolded spectrum is shown on the right panel of Fig. 1 (we considered about 10 5 levels in the bulk of the spectrum). We observe that P(x) is maximum at x = 0 indicating level clustering, and it is likely to be the negative exponential P(x) e −x characteristic of the Poisson process. We have also studied other spin models obtaining similar results. This was to be expected since the XY model and its variants are integrable. Poisson statistics have also been numerically observed in previous works for other spin systems integrable by Bethe ansatz, including the Heisenberg chain, the t-J model, and the Hubbard model. See, e.g., Ref. 29.
We have performed the same investigation for random quadratic forms with independent Gaussian coefficients (see Sec. VI), where the elementary excitations λ k,n of the normal modes are distributed as the singular values of the real Ginibre ensemble (61). Our findings are reported in Fig. 2. Again, the Gaussian limit (60) is a convincing approximation of the numerical level density even for moderate values of n (left panel). The level spacing in the unfolded spectrum (about 10 5 levels in the bulk) is well described by a negative exponential. Note that the elementary excitations λ k,n repel as the eigenvalues of random matrices (see Eq. (62)); nevertheless, the energy levels E k ,n are given by the subset sums of the λ k,n 's, and this structure dominates the repulsion and enhances the presence of small gaps. At first, this result may be surprising for those working in the field of random matrices or spectral theory of disordered systems. For generic chaotic systems, one usually expects level repulsion. We felt natural to provide a theoretical argument to explain the "lack of repulsion" for disordered quasifree fermions.
As argued theoretically by Berry and Tabor, 8 the energy spectrum of a classically integrable Hamiltonian system represents a sequence of completely uncorrelated numbers and the spectral fluctuations obey Poissonian statistics. The original argument in Ref. 8 is based on the fact that for integrable systems, it is possible to perform a canonical transformation into action-angle coordinates. The semiclassical approximation consists in quantizing the action variables so that the quantum energy levels of a classically integrable system are given by the classical Hamiltonian evaluated at points of a lattice (in some cases this quantization rule is exact). Therefore, the level spacings or, more generally, the number statistics of energy levels are related to the problem of counting the number of lattice points enclosed by the Hamiltonian level sets. A computation based on the Poisson summation formula then suggests that P(x) exp(−x) for generic integrable systems. This scheme applies only to "generic" systems, and some notable exceptions are quite well known.
Later, this way of reasoning has been extended beyond Hamiltonian mechanics. For instance, the standard argument for Poisson statistics in the case of spin integrable models is as follows. 29 If a Bethe ansatz holds, the energy levels of the systems are characterised by a set of quasimomenta (that reduce to real momenta for noninteracting spin systems). Typically, these quasi-momenta are the solutions of a set of non-linear equations, and therefore, the possible quasimomenta are likely to repel one another, namely, they lie on a quasilattice. The level statistics again reduces to the statistics of the lattice positions, and the same argument as Ref. 8 leads to Poisson statistics.