text:
# Spectral ergodicity and normal modes in ensembles of sparse matrices

## I Introduction

Random matrix theory (RMT) has shown itself to be useful in modeling a wide variety of physical systems . Sparse random matrices, which are characterized by a considerable fraction of vanishing or negligible small matrix elements, have attracted particular interest in recent years. Applications of sparse matrices can be found, for instance, in disordered systems, quantum chaotic systems or many-body systems . There are few analytical investigations of sparse matrices, and these are often restricted to the subclass of random band matrices, i.e. matrices where the off-diagonal elements vanish for a sufficiently large distance from the main diagonal. More general classes of sparse matrices must be treated mainly by numerical methods.

In RMT, spectral quantities are usually calculated as averages over the matrix ensemble. Experimental data, however, normally consists of only a single spectrum, and the corresponding quantities must therefore be calculated as a running average over the energy. Comparison of experiment and theory is facilitated by the spectral ergodicity hypothesis for matrix ensembles: An average over an ensemble of random matrices provides the same result as an average over energies for a single element of this ensemble in the limit where the dimension of the matrices is large. In the case of full (i.e. non-sparse) random matrices only, spectral ergodicity can be proved rigorously . The question of whether a given matrix ensemble obeys spectral ergodicity is not merely academic; it has a direct bearing on the ability to compare experimental results with theoretical expectations.

In the context of disordered mesoscopic systems, the issue of ergodicity was first raised in Ref. . Here the ergodicity hypothesis states that averages over all realizations of the disorder potential yield the same values for an observable as averages over an interval variable such as energy or an applied magnetic field. For such systems – we think of a single electron in a crystal with random impurities – the Hamiltonian is modeled by a sparse random matrix via a spatial discretization on a lattice. The disorder potential is then represented by an appropriate choice of the matrix elements as random variables. The diagonal elements represent the on-site energies; the off-diagonal elements represent the coupling between the sites. Assuming that only neighboring sites are coupled, one obtains a Hamiltonian of sparse matrix form. One-dimensional systems, e.g. long and thin wires, are described by random band matrices . In higher dimensions, the inevitable presence of side bands renders analytical treatment difficult.

The question of the ergodicity of mesoscopic disordered systems is important for the extraction of the Thouless energy from spectral data. The Thouless energy, $`E_\mathrm{c}`$, is an intrinsic energy scale in disordered systems  which is essentially the inverse of the time required for a particle to move diffusively through the sample as a consequence of multiple scatterings by impurities. The Thouless energy is measured in terms of the mean level spacing, $`\mathrm{\Delta }`$, and can be extracted from the fluctuations in the eigenvalue spectrum of the corresponding Hamiltonian. Here the mean level spacing $`\mathrm{\Delta }`$ is the inverse of the spectral density, $`\rho `$, and is a function of energy. Furthermore, the Thouless energy (in units of the mean level spacing) is the dimensionless conductance, $`g=E_\mathrm{c}/\mathrm{\Delta }`$, a parameter which is sufficient to describe important qualitative features of disordered mesoscopic systems. For $`g1`$ the system is a good conductor, and the eigenfunctions are delocalized. Similarly, $`g1`$ characterizes an insulator with localized eigenfunctions. A transition between both regimes is observed for $`g1`$. Moreover, spectral fluctuations for level spacings smaller than $`E_\mathrm{c}/\mathrm{\Delta }`$ obey Wigner-Dyson statistics while for larger level spacings corrections to RMT behavior arise. The actual value of $`g`$ depends on the degree of disorder and on the spatial dimension of the system. In theoretical investigations, $`g`$ is usually calculated from an ensemble average. In experiments, however, it must frequently be obtained from a single sample, i.e. a single realization of the disorder potential. The ergodicity hypothesis requires that both values of $`g`$ are the same. Due to the observation that the size of sample-to-sample fluctuations is of the same order as the fluctuations induced by changes in the energy or magnetic field , this hypothesis appears to be valid in such systems. Since any particular sample has a definite conductance $`g`$, which can be measured independent of the existence of other samples, the ergodicity hypothesis in disordered systems is tantamount to the statement that the Thouless energy can be extracted from a single member of the ensemble.

An apparent counter-example to the spectral ergodicity hypothesis has been found in another class of systems. Complex many-body systems such as atomic nuclei can be modeled in the framework of random-matrix theory by the so called two-body random ensemble (TBRE) . In this model, fermionic particles are coupled by a stochastic two-body interaction and distributed over a set of single-particle states. The many-body Hamiltonian is then represented by a matrix in the basis of the many-body states constructed as Slater determinants of the single-particle wave functions. Due to conservation of total angular momentum, $`J`$, and total isospin, $`T`$, the matrix has block–diagonal structure. More importantly, the two–body nature of interaction ensures that each block with given $`J`$ and $`T`$ is sparse. Numerical results indicate that the TBRE is non-ergodic in the sense that spectral and ensemble averages do not coincide . One finds the results of spectral averaging in agreement with experiment .

A second counter-example to the spectral ergodicity hypothesis can be found in the spectral fluctuations of the lattice QCD Dirac operator . Again, recent developments suggest a connection to disordered systems. Using formal analogies, it has been argued that there exists an equivalent of the Thouless energy in QCD . This energy scale can be associated with properties of the QCD vacuum related to the spontaneous breaking of chiral symmetry. This resembles the situation in solid states physics where $`E_\mathrm{c}`$ is related to the ground state properties of the solid, i.e. $`g`$. These expectations have been confirmed by numerical lattice gauge simulations . Moreover, generalization of the analytic treatments from solid state physics to include chiral symmetry permits a more rigorous derivation of the analogies between disordered systems and QCD . The essential result is that the spectral fluctuations of the QCD Dirac operator exhibit RMT behavior up to a certain scale, the equivalent of a Thouless energy. Beyond this scale, corrections to the universal statistics arise. The spectral data that is used to extract the Thouless energy from Dirac spectra is obtained from an ensemble average over gauge field configurations. This is the most natural approach in QCD, since the QCD partition function is defined as a path integral, i.e. an ensemble average. If spectral averaging is performed nevertheless, an unexpected result is obtained. The energy scale seen in the ensemble average of the spectral statistics disappears, and RMT fluctuations are observed to almost arbitrarily large scales . This observed violation of spectral ergodicity may not be disturbing given the intrinsic nature of ensemble averaging in QCD. It is currently impossible to say whether this represents a fundamental difference between QCD and disordered systems or whether it may be possible to find similar effects in disordered systems.

We have mentioned three examples of physical systems in which sparse matrix models apply. In disordered mesoscopic systems the ergodicity hypothesis is an essential and apparently valid assumption. The two counter-examples indicate that spectral ergodicity is not necessarily respected. In each of these cases, however, there are clear physical arguments indicating which averaging procedure is appropriate. For the TBRE, the spectral average is natural since one deals with spectra obtained from single nuclei. In lattice QCD, on the other hand, the ensemble average is intrinsic to the system from its very construction. Nevertheless, it seems important to determine which spectral properties of these systems are responsible for the difference between spectral and ensemble averaging and to understand whether similar effects, disproving the ergodicity hypothesis, can also be found in disordered systems. One can also imagine new classes of random matrix models which violate spectral ergodicity but for which physical arguments in favor of a specific averaging procedure are not readily available. A better understanding of the mechanisms responsible for the breakdown of ergodicity could prove to be essential in such systems. The three examples considered here are, at first sight, extremely different in character. As noted, their most obvious common feature is the sparsity of the matrix representations of their quantum mechanical Hamiltonian operators. Therefore, a detailed investigation of the influence of sparsity on matrix models and on the possible breaking of ergodicity in them seems to be desirable.

Investigations of the TBRE have suggested that the differences between spectral and ensemble averages are related to fluctuations of the average spectral density over the ensemble . In this light, the problem is how to determine the average spectral density. There is little ambiguity when constructing the ensemble average. By averaging over a sufficiently large number of members of the ensemble, it is easy to determine the spectral density in an energy interval of any given size. The only delicacy in this procedure lies in the size of the energy intervals chosen. We are interested in constructing an average local spectral density in the thermodynamic limit of large matrices. In some cases, such as the Gaussian ensembles, this double limit of vanishing interval and infinite particle number is unproblematic. This is not necessarily the case. Concerning spectral averaging, it is necessary to adopt some kind of local smoothing. Unfortunately, there is no rigorous and well defined procedure to accomplish this task, which is usually regarded as an annoying and cumbersome technical detail. A primary aim of this paper is, however, to demonstrate that it is precisely the definition of smoothing which provides the key to understanding the apparent lack of ergodicity in some systems.

In considering the effect of fluctuations on averaging procedures, it is essential to recognize that the fluctuations in individual eigenvalues of a random matrix are not statistically independent. Indeed, it can be shown that the statistically independent fluctuations involve the collective motion of literally all eigenvalues in the spectrum. The statistically independent normal modes of the spectrum provide a suitable tool for describing this collective motion . As we will show, the nature of these normal modes is extremely simple, and they can be regarded as plane waves with a well defined wave length. These normal modes can be crudely divided into two classes. Short wave fluctuations should provide information regarding universal spectral properties, and an appropriate averaging procedure should retain their effects. On the other hand, long wave length modes describe model-dependent, non-universal physics and should be eliminated by the averaging procedure. The art of averaging thus lies in establishing a physically sensible division between short and long wave lengths. (In this sense, the challenge of constructing an appropriate spectral average is quite similar to the task of defining an appropriate energy interval when performing ensemble averages.) Once this division has been made, the question of the validity of the ergodic hypothesis can be answered readily by considering the mean square amplitude of the normal modes as a function of wave length. In some cases, these mean square amplitudes have a simple functional dependence on the wavelength, which applies with equal validity to fluctuations of all wave lengths. This is the case for the familiar Gaussian ensembles of non-sparse random matrices, where mean square amplitudes are linearly proportional to the wave length. In such cases, the details of the spectral averaging process are irrelevant, and random matrix ergodicity is respected. In other cases, including the case of sparse random matrices, the situation is more delicate. There we will encounter a qualitative difference which, not surprisingly, suggests a clear distinction between the “softness” of long wave length fluctuations and the relative “rigidity” of short wave length modes. Spectral averaging as commonly employed tends to be more efficient than ensemble averaging in eliminating these soft long wave length fluctuations. The result of this argument will provide a simple and natural explanation of the differences between ensemble and spectral averaging found in these systems. Further, it will enable us to understand why the results of spectral averaging show greater consistency with the familiar results of the non-sparse Gaussian ensembles. We shall finally see that the apparent breakdown of random matrix ergodicity is, in fact, a false puzzle. If we ask the same physical question, the two averaging methods will provide us with the same answer. If, as can be the case for sparse random matrices, ensemble and spectral averaging probe different and conflicting aspects of the system, we should not be startled to obtain different results.

The paper is organized as follows. We introduce our matrix model in Sec. II. The relevant parameters of the model are discussed in Sec. III, which is followed by a brief discussion of the spectral density as a function of these parameters in Sec. IV. In Sec. V, we present a detailed analysis of the spectral correlations obtained from the ensemble average. This includes short range as well as long range correlations. A critical scale is found in the spectral statistics, which is interpreted in terms of a Thouless energy. The dependence of this scale on the parameters of the matrix model is discussed. The findings from the ensemble average are contrasted with those obtained from spectral averaging in Sec. VI. The spectral observables obtained from the two averaging procedures do not agree. Moreover, the results from spectral averaging depend on how the local smoothing of the spectral density is performed. This difference is qualitatively explained by the collective motion of the eigenvalues. Normal modes are thus introduced in Sec. VII to provide a suitable mathematical description of this collective motion. They will provide us with a natural explanation of the differences between the results of Sects. V and VI and help us to understand that there is no genuine violation of random matrix ergodicity. We will offer a summary and conclusions in Sec. VIII.

## II Sparse random matrix model

We are interested in some generic features of systems such as disordered systems, TBRE, and lattice QCD, which can all be described with sparse random matrices. Despite the physical differences in these systems, their spectral properties show similar features. We will therefore concentrate on the influence of sparsity on the spectral properties. We shall not attempt to model the fine structure of the TBRE or the multi-band structure of a $`d`$-dimensional disordered system. We will not include the combinatorial correlations of the TBRE. We will certainly not attempt to model the complicated structure of non-abelian lattice gauge theories.

Our matrix model consists of an ensemble of sparse real symmetric matrices. We introduce the sparsity $`\alpha `$, which is the fraction of the $`N(N1)/2`$ independent off-diagonal matrix elements, chosen to be non-vanishing. All diagonal elements are kept non-zero. The non-vanishing matrix elements $`H_{ij}`$ are chosen independently and at random from a Gaussian distribution with mean $`0`$ and a variance which is for the diagonal elements $`\sqrt{2}\sigma `$ and for the off-diagonal $`\sigma `$,

$$P(H_{ij})=\frac{1}{\sqrt{2\pi \sigma _{ij}^2}}\mathrm{exp}\left(\frac{H_{ij}^2}{2\sigma _{ij}^2}\right),$$
(1)
with $`\sigma _{ij}^2=1+\delta _{ij}`$. This ensures that GOE behavior is recovered in the limiting case $`\alpha =1`$.

Since we are not concerned with the effects of any particular configuration of the non-vanishing off-diagonal elements, we choose their positions at random and independently for each matrix of the ensemble. We have, however, verified that our results do not change if one instead maintains randomly chosen but fixed positions for the non-vanishing elements.

## III Effective dimension

Our matrix model depends only on two parameters, the matrix size $`N`$ and the sparsity $`\alpha `$. While analytic treatment is difficult for arbitrary $`\alpha `$, it is possible to find a qualitative description of the dependence of spectral correlations on the parameters of the model. In particular, comparison with the properties of disordered systems suggests that there is only one relevant parameter, the effective dimension.

Consider a particle moving in a $`d`$-dimensional disordered medium . Suppose the strength of the disorder is chosen so that the motion of the particle is diffusive. It follows from the diffusion equation appropriate for the system that the mean square distance, $`x^2`$, traversed by the particle in time $`t`$ is proportional to $`t`$

$$x^2=2d𝒟t,$$
(2)
where $`𝒟`$ is the diffusion constant. If the system is restricted to a finite volume $`V`$, the particle will have explored the entire system for times larger than the diffusion time. For such times the probability of finding the particle at any position $`x`$ is everywhere equal. The diffusion time is given by

$$t_\mathrm{d}=\frac{V^{2/d}}{2d𝒟},$$
(3)
which follows from Eq. (2) if the mean square distance equals the linear size of the system, $`\sqrt{x^2}=V^{1/d}`$. The Thouless energy $`E_\mathrm{c}`$ is simply the inverse of the diffusion time $`E_\mathrm{c}1/t_\mathrm{d}`$. This classical quantity can be observed in the correlations of the eigenvalues of the quantum mechanical spectrum. Spectral correlations are commonly examined on the scale of the mean level spacing $`\mathrm{\Delta }`$, which scales like $`\mathrm{\Delta }1/V`$. Combining the above arguments yields a dimensionless critical energy, $`L_\mathrm{c}`$, given by

$$L_\mathrm{c}=E_\mathrm{c}/\mathrm{\Delta }d𝒟V^{12/d}.$$
(4)
If the levels have a separation less than $`L_\mathrm{c}`$, their correlations follow the predictions of random matrix theory. If their separation is greater than $`L_\mathrm{c}`$, corrections to these universal fluctuations arise.

For $`d<2`$, states are localized for any strength of the disorder . The arguments above cannot be applied since the motion of the particle is not described by a diffusion equation, which is an essential assumption in arriving at Eq. (2). For $`d>2`$, states will be delocalized for sufficiently small disorder. In that case Eq. (4) is valid, and $`L_\mathrm{c}`$ can be determined from the spectral statistics.

In order to establish a link to our matrix model, we think of a suitable discretized lattice version of the corresponding Hamiltonian, $`H`$. In the diffusive regime, the above arguments are also valid for a lattice realization of the Hamiltonian. One now expects to find a critical energy, $`L_\mathrm{c}`$, which scales with the volume (now given by the matrix dimension $`N`$) in a manner similar to Eq. (4). Thus, we expect that $`L_\mathrm{c}(\alpha ,N)𝒞(\alpha ,N)N^{12/d(\alpha ,N)}`$ for some $`𝒞(\alpha ,N)`$ and $`d(\alpha ,N)`$.

For a given $`N`$ and $`d`$, it is easy to estimate $`\alpha `$. There are $`N(N1)/2`$ independent off-diagonal elements. For each lattice there are $`d`$ couplings to neighboring sites in one direction, one for each dimension. Thus, there are $`dN`$ non-vanishing matrix elements out of the total $`N(N1)/2`$. The sparsity is then

$$\alpha (N,d)=\frac{dN}{N(N1)/2}=\frac{2d}{N1}.$$
(5)
Solving for $`d`$, one obtains an effective dimension $`d(\alpha ,N)`$ for any $`\alpha `$ and $`N`$.

However, the argument leading to Eq. (5) is qualitative, and several reservations should be mentioned. Firstly, the coupling of sites need not be restricted to nearest neighbors, and this can lead to different factors dependent on the number of couplings. Secondly, the above arguments are strictly valid only for cubic lattices; other lattice geometries could lead to additional factors. It remains nevertheless that the number of couplings grows like $`dN`$ and the number of matrix elements grows like $`N^2`$ for large $`N`$. Thus, $`\alpha `$ should be proportional to $`d/N`$. We therefore introduce an effective dimension defined as

$$d_{\mathrm{eff}}=\alpha N$$
(6)
without additional constants. By construction, $`0\alpha 1`$ and thus $`0d_{\mathrm{eff}}N`$. The effective dimension can also be understood as twice the average number of independent off-diagonal elements per row or per column. For sufficiently large dimensions, spectral statistics should possess a certain scale $`L_\mathrm{c}`$ below which they follow RMT predictions and above which stronger non-RMT fluctuations appear. Thus, we expect that $`L_\mathrm{c}`$ should have the form

$$L_\mathrm{c}(N,d_{\mathrm{eff}})=\{\begin{array}{cc}c(d_{\mathrm{eff}})N^{\eta (d_{\mathrm{eff}})}\hfill & d_{\mathrm{eff}}>d_\mathrm{c}\hfill \\ & \\ 0\hfill & d_{\mathrm{eff}}<d_\mathrm{c}.\hfill \end{array}$$
(7)
For $`d_{\mathrm{eff}}<d_\mathrm{c}`$, the states should be localized. In this region, $`L_\mathrm{c}1`$, and the spectral statistics should be those of a Poisson distribution. This is obviously true in the limit $`d_{\mathrm{eff}}0`$, where the coupling between diagonal elements vanishes. The eigenvalues of $`H`$ are then the randomly distributed diagonal elements, which obey the Poisson spectral statistics of an uncorrelated sequence of levels. In the opposite limit, $`d_{\mathrm{eff}}>d_\mathrm{c}`$, states are expected to be delocalized. Spectral statistics will have the Wigner-Dyson form up to $`L_\mathrm{c}`$. In the extreme limit $`d_{\mathrm{eff}}=N`$, $`L_\mathrm{c}`$ is equal to $`N`$ since this limit corresponds to the original GOE ensemble. At a certain value there has to be necessarily a transition between the both limits. The exponent should obey

$$\eta (d_{\mathrm{eff}})=ab/d_{\mathrm{eff}},$$
(8)
where the constants $`a`$ and $`b`$ are independent of $`d_{\mathrm{eff}}`$, $`\alpha `$ and $`N`$. Since experience does not indicate that spectral properties are violently sensitive to the sparsity, we expect that these parameters will be related so that $`\eta (d_\mathrm{c})0`$. Thus, we expect that

$$d_\mathrm{c}b/a.$$
(9)
The constant $`c(d_{\mathrm{eff}})`$ in Eq. (7) should depend only on the effective dimensionality of the system in the case of fixed disorder strength. This can be understood from the definition of the diffusion constant as the second moment of the probability density  in the infinite volume limit, which depends on the dimension and the degree of disorder. The latter is kept fixed in our model.

## IV Spectral density

For arbitrary sparsity, an analytic description of the spectral density of the matrix model is difficult . In the limiting cases $`d_{\mathrm{eff}}=0`$ and $`d_{\mathrm{eff}}=N`$, elementary analytic expressions are readily available.

For $`d_{\mathrm{eff}}=0`$, the Hamiltonian matrix reduces to a diagonal matrix with independent, Gaussian distributed eigenvalues. From Eq. (1) one immediately obtains

$$\underset{d_{\mathrm{eff}}0}{lim}\rho (E)=\frac{N}{2\sqrt{\pi }}\mathrm{exp}\left(\frac{E^2}{4}\right).$$
(10)
Empirical corrections to the Gaussian shape for small but non-zero $`d_{\mathrm{eff}}`$ were calculated in Ref. . On the other hand, the matrix model for $`d_{\mathrm{eff}}=N`$ coincides with the standard Gaussian orthogonal ensemble. The spectral density is then known to have a semi-circular shape ,

$$\underset{d_{\mathrm{eff}}N}{lim}\rho (E)=\frac{1}{2\pi }\sqrt{4NE^2}.$$
(11)
The spectral density interpolates smoothly between these limits as $`d_{\mathrm{eff}}`$ is varied. The spectral density is conventionally normalized to the total number of eigenvalues, $`_{\mathrm{}}^{\mathrm{}}𝑑E\rho (E)=N`$. Fig. 1 shows $`\rho (E)/N`$ for fixed $`d_{\mathrm{eff}}`$ and varying $`N`$. The ratio $`\rho (E)/N`$ depends solely on $`d_{\mathrm{eff}}`$. The overall shape evolves smoothly from an approximate semi-circle to an approximate Gaussian. This is consistent with the findings from many-particle spectra .

Knowledge of the spectral density is crucial in the statistical analysis of spectra because spectral correlations must be investigated on a scale set by the local mean level spacing, $`\mathrm{\Delta }`$, which is given by $`\mathrm{\Delta }(E)=1/\rho (E)`$. The transformation

$$\xi =_{\mathrm{}}^E𝑑E^{}\rho (E^{})$$
(12)
leads to a dimensionless energy variable, $`\xi `$, which is the energy measured in terms of the mean level spacing. This procedure is referred to as “unfolding”. The purpose of this unfolding is to eliminate model-dependent macroscopic variations in the spectral density in order to reveal underlying universal spectral correlations.

The spectral density can be determined in two ways which differ in principle. The first is to calculate $`\mathrm{\Delta }`$ from an ensemble average. The second is to smooth the actual spectral density by local averaging. The latter approach is the conventional way of determining $`\mathrm{\Delta }`$ and has been applied to a variety of systems . We will refer to this method as self-unfolding or spectral unfolding. Obviously, when we have only a limited number of experimental samples at our disposal, spectral unfolding is the only approach available. While the spectral densities obtained with these methods are superficially similar, this can be misleading. In TBRE  and lattice QCD , these two spectral densities lead to very different results. Ensemble unfolded spectra show a critical scale beyond which spectral fluctuations are no longer described by RMT. By contrast, self-unfolded spectra show fluctuations which are in excellent agreement with RMT expectations on all energy scales.

It is important to emphasize that these observed differences are due to the way in which the spectra are unfolded and not to the way in which spectral correlations within the ensemble of unfolded spectra are calculated. This fact is frequently ignored and contradicts the common view that unfolding is a purely technical procedure with no physical content. We will discuss this point at some length.

In the following, we denote ensemble averages by a bar, $`\overline{(\mathrm{})}`$. Spectral averaged quantities will be denoted by $`\mathrm{}`$.

## V Spectral observables from ensemble unfolding

In this Section, we investigate several spectral observables using ensemble unfolding. We will concentrate on the dependence of these observables on the effective dimension, $`d_{\mathrm{eff}}`$. In Sec. V A, we study the transition from Poisson to Wigner statistics exhibited by our matrix model in terms of the nearest neighbor spacing distribution. Long range spectral correlations and the functional dependence of the Thouless energy on $`N`$ and $`d_{\mathrm{eff}}`$ are investigated in Sec. V B.

### A Short range correlations

The nearest neighbor spacing distribution, $`P(s)`$, is the probability density of finding two adjacent levels at a distance $`s=\xi _{i+1}\xi _i`$. For a sequence of uncorrelated levels, i.e. the Poisson case, this distribution is simply $`P_\mathrm{P}(s)=\mathrm{exp}(s)`$. In the case of random matrices, the nearest neighbor distribution is well approximated by the Wigner surmise, which reads $`P_{\mathrm{WD}}(s)=\pi s\mathrm{exp}(\pi s^2/4)/2`$ for real symmetric matrices.

Fig. 2 shows the nearest neighbor distribution for fixed $`N`$ and varying $`d_{\mathrm{eff}}`$. When $`d_{\mathrm{eff}}>4`$, the data is described well by the Wigner surmise. A transition between Wigner and Poisson forms is seen in the vicinity of $`d_{\mathrm{eff}}=2`$. In the limit of low effective dimension, $`P(s)`$ is accurately described by Poisson statistics. This is in agreement with theoretical considerations , which suggest that there should be a transition when the average number of independent non-vanishing off-diagonal elements is of order one per row.

In order to obtain a more quantitative analysis of the transition of $`P(s)`$ observed in Fig. 2, we evaluate the integral of the tail of $`P(s)`$, $`A=_{s_0}^{\mathrm{}}P(s)𝑑s`$ , which is free from binning effects, for various values of $`d_{\mathrm{eff}}`$. Here $`s_02.002`$ is the value at which the curves of Poisson and Wigner form cross. We the Wigner-Dyson and Poisson values as $`A_{\mathrm{WD}}`$ and $`A_\mathrm{P}`$, respectively. The parameter

$$\gamma (d_{\mathrm{eff}})=\frac{A(d_{\mathrm{eff}})A_{\mathrm{WD}}}{A_\mathrm{P}A_{\mathrm{WD}}}$$
(13)
then gives us a quantitative description of the statistics of the spectra with values $`0\gamma 1`$. The transition is smooth for all finite matrix sizes but becomes sharper as $`N`$ increases. A comparison of $`\gamma (d_{\mathrm{eff}})`$ for different sizes of the system reveals that all curves cross at the same value for the effective dimension where the size effect changes its sign .

In Fig. 3, we plot $`\gamma (d_{\mathrm{eff}})`$ as a function of the effective dimension for three different values of $`N`$. The transition from Poisson to Wigner–Dyson spectral fluctuations is evident. The three curves saturate at the expected values $`\gamma =1`$ for small values of $`d_{\mathrm{eff}}`$ and $`\gamma =0`$ at large $`d_{\mathrm{eff}}`$. As expected, the transition becomes sharper as $`N`$ increases. It is also clear from the figure that all curves cross at $`d_{\mathrm{eff}}2`$. This suggests the existence of a critical dimension for the transition between Poisson and Wigner type fluctuations at $`d_{\mathrm{eff}}2`$.

### B Long range correlations

The correlator $`P(s)`$ and the related parameter $`\gamma `$ indicate that there is a qualitative change in the short-range spectral fluctuations in the vicinity of $`d_{\mathrm{eff}}=2`$. In this Section, we will extend this analysis to long-range fluctuations, i.e. to scales which are significantly larger than the mean level spacing. We do this by studying the number variance. All results in this Section have been obtained from the fluctuations of the ensemble unfolded eigenvalues. A similar analysis of the self-unfolded eigenvalues and the significant changes in the properties of long-range fluctuations found there will be discussed in Sec. VI.

The number variance measures fluctuations in the total number of levels found in an energy interval $`[L_0L/2,L_0+L/2]`$. In this Section $`L=_{\mathrm{}}^E𝑑E^{}\overline{\rho }(E^{})`$ is the ensemble unfolded energy. The definition of the unfolding process ensures that the average number of levels in this interval is $`L`$ independent of $`L_0`$. The number variance is then given by

$$\overline{\mathrm{\Sigma }}_{L_0}^2(L)=\overline{(n_{L_0}(L)\overline{n}_{L_0}(L))^2}.$$
(14)
The most conservative way to calculate observables such as the number variance is to perform the ensemble average with fixed energy, $`L_0`$. We have, however, verified that our results are independent of $`L_0`$ provided that we avoid the edges of the spectrum. This property, familiar from the Gaussian random matrix ensembles, is called translational invariance. Numerical investigation shows that translational invariance is quantitatively reliable except for some $`10\%`$ of the total number of eigenvalues in the vicinity of each of the edges. We will restrict our attention to the translationally invariant bulk of the spectra. Given translational invariance, we will omit the subscript $`L_0`$ in the following.

A sequence of uncorrelated levels, the Poisson case, gives rise to a strictly linear number variance, $`\mathrm{\Sigma }^2(L)_{\mathrm{Poisson}}=L`$. By contrast, random matrices have much stronger correlations, and the number variance grows only logarithmically for large $`L`$. The asymptotic form for the GOE result for $`L1`$ is given by $`\mathrm{\Sigma }^2(L)_{\mathrm{GOE}}=2/\pi ^2[\mathrm{ln}(2\pi L)+\gamma +1\pi ^2/8]`$, where $`\gamma `$ is Euler’s constant .

Fig. 4 shows the number variance as a function of $`d_{\mathrm{eff}}`$. The matrix dimension is fixed at $`N=1000`$. The number variance is calculated from the central interval $`[L/2,L/2]`$ averaged over the ensemble. No spectral averaging is performed. Note that the same ensemble of matrices was used for the calculation of the number variance at each value of $`L`$. This inevitably results in strongly correlated statistical uncertainties, which are clearly visible in the figure. Above the critical dimension, $`d_\mathrm{c}2`$, the number variance is described by GOE statistics up to a certain scale $`L_\mathrm{c}(N,d_{\mathrm{eff}})`$. For $`L>L_\mathrm{c}(N,d_{\mathrm{eff}})`$, fluctuations become stronger and can no longer be described by the GOE. This can be seen in the upper row of Fig. 4. The critical scale $`L_\mathrm{c}(N,d_{\mathrm{eff}})`$ decreases as $`d_{\mathrm{eff}}`$ approaches $`d_\mathrm{c}`$ from above. The change in the spectral statistics as one crosses $`d_\mathrm{c}`$ is clearly visible in the lower row of Fig. 4. Well below $`d_\mathrm{c}`$, fluctuations show the purely linear behavior of the Poisson distribution. This is similar to the transition seen in the nearest neighbor spacing distribution of Fig. 2. The critical scale $`L_\mathrm{c}(N,d_{\mathrm{eff}})`$ is of order one in the vicinity of $`d_{\mathrm{eff}}2`$ and gives rise to the intermediate statistics seen in the nearest neighbor spacing distribution.

The scale $`L_\mathrm{c}(N,d_{\mathrm{eff}})`$ is not uniquely defined. We elected to use the following definition:

$$\frac{|\mathrm{\Sigma }^2(L,N,d_{\mathrm{eff}})_{\mathrm{data}}\mathrm{\Sigma }^2(L)_{\mathrm{GOE}}|}{\mathrm{\Sigma }^2(L)_{\mathrm{GOE}}}>\epsilon \mathrm{for}L>L_\mathrm{c}^{(\epsilon )}(N,d_{\mathrm{eff}}).$$
(15)
The precise value of $`L_\mathrm{c}^{(\epsilon )}(N,d_{\mathrm{eff}})`$ depends on the choice of $`\epsilon `$. Moreover, this definition leads to a non-zero value of $`L_\mathrm{c}^{(\epsilon )}(N,d_{\mathrm{eff}})`$ whenever $`\epsilon >0`$ — even when the data obeys Poisson statistics. The minimal value of $`L_\mathrm{c}^{(\epsilon )}(N,0)`$ in that case follows from the solution of Eq. (15) with $`\mathrm{\Sigma }^2(L,N,d_{\mathrm{eff}})_{\mathrm{data}}=L`$, from which follows that

$$\frac{L_\mathrm{c}^{(\epsilon )}(N,0)}{\mathrm{\Sigma }^2(L_\mathrm{c}^{(\epsilon )}(N,0))_{\mathrm{GOE}}}=1+\epsilon .$$
(16)
Using the analytic form of the GOE number variance for small $`\epsilon `$, it is easy to show that the solution to Eq. (16) for a pure Poisson distribution is $`L_\mathrm{c}^{(\epsilon )}(N,0)=\epsilon `$ for small $`\epsilon `$. Thus, the observation of $`L_\mathrm{c}^{(\epsilon )}(N,d_{\mathrm{eff}})\epsilon `$ implies that $`L_\mathrm{c}(N,d_{\mathrm{eff}})=0`$.

The qualitative arguments in Sec. III suggest that $`L_\mathrm{c}`$ should have a power-law dependence on $`N`$ for $`d_{\mathrm{eff}}>d_\mathrm{c}`$. Both the exponent and the effective diffusion coefficient should depend only on $`d_{\mathrm{eff}}`$. Thus, we calculate $`L_\mathrm{c}^{(\epsilon )}(N,d_{\mathrm{eff}})`$ for various $`N`$ and fixed $`d_{\mathrm{eff}}`$. Some examples are shown in Fig. 5. They indicate a clear power-law dependence on $`N`$. The value of $`L_\mathrm{c}^{(\epsilon )}(N,d_{\mathrm{eff}})`$ is seen to decrease with decreasing $`d_{\mathrm{eff}}`$.

The slope in Fig. 5 is simply the scaling exponent, $`\eta (d_{\mathrm{eff}})`$, introduced in Eqs. (7) and (8). It can be extracted from the data by a linear fit. The scaling exponent as a function of the effective dimension is shown in Fig. 6 as a function of $`\epsilon `$. The error bars represent the variance of the linear fit. The exponent $`\eta (d_{\mathrm{eff}})`$ drops to zero at about $`d_{\mathrm{eff}}=2`$ in agreement with our expectations. For $`d_{\mathrm{eff}}<2`$, the spectral statistics are of Poisson form, and $`L_\mathrm{c}^{(\epsilon )}(N,d_{\mathrm{eff}})\epsilon `$ independent of $`N`$ as suggested by Eq. (16). In contrast to $`\eta (d_{\mathrm{eff}})`$, the diffusion coefficient in Eq. (7), now denoted as $`c^{(\epsilon )}(d_{\mathrm{eff}})`$, does depend on $`\epsilon `$. In the Poisson regime, the original definition of the diffusion coefficient requires $`c(0)=0`$. However, given the construction Eq. (15), the diffusion coefficient is non-zero even for $`d_{\mathrm{eff}}<d_\mathrm{c}`$. In the limit $`d_{\mathrm{eff}}=0`$, we have $`\eta (0)=0`$ and therefore $`c^{(\epsilon )}(0)=L_\mathrm{c}^{(\epsilon )}(N,0)`$. We thus normalize the diffusion coefficient obtained from the linear fit with $`c^{(\epsilon )}(0)`$. The numerical results are shown in Fig. 7. The normalized diffusion coefficient is found to be remarkably insensitive to $`\epsilon `$. Thus, we conclude from our matrix model that

$$\frac{L_\mathrm{c}^{(\epsilon )}(N,d_{\mathrm{eff}})}{L_\mathrm{c}^{(\epsilon )}(N,0)}=C(d_{\mathrm{eff}})N^{\eta (d_{\mathrm{eff}})}d_{\mathrm{eff}}>d_\mathrm{c}.$$
(17)
The right hand side of this equation is independent of $`\epsilon `$, which only influences the overall scale. The scaling exponent $`\eta (d_{\mathrm{eff}})`$ and the normalized diffusion coefficient $`C(d_{\mathrm{eff}})=c^{(\epsilon )}(d_{\mathrm{eff}})/c^{(\epsilon )}(0)`$ depend only on the effective dimension. Below the critical dimension we have $`\eta (d_{\mathrm{eff}}<d_\mathrm{c})=0`$. Above the critical dimension, the suggested functional dependence, Eq. (8), is $`\eta (d_{\mathrm{eff}}>d_\mathrm{c})=ab/d_{\mathrm{eff}}`$. From the analysis of the short range correlators, we have $`d_\mathrm{c}=2`$, which suggests that $`a=b/2`$. Performing a one parameter fit to the data for $`d_{\mathrm{eff}}>2`$, we actually find $`a=0.51\pm 0.05`$. The function $`\eta (d_{\mathrm{eff}})=1/21/d_{\mathrm{eff}}`$ is also shown in Fig. 6 and is seen to be in reasonable agreement with the data.

## VI Spectral observables from self-unfolding

We now perform an analysis similar to that of the previous section in which ensemble unfolding is replaced by self-unfolding $`L=_{\mathrm{}}^E𝑑E^{}\rho (E^{})`$. Here, $`\rho (E)`$ is the spectral density obtained from the local smoothing of each member of the ensemble. The number variance is, in analogy to Eq. (14), constructed as

$$\mathrm{\Sigma }^2(L)=\mathrm{\Sigma }_{L_0}^2(L)=(n_{L_0}n_{L_0})^2.$$
(18)
The average is performed over various intervals $`[L_0L/2,L_0+L/2]`$ with fixed interval length, in one spectrum. The final number variance is then obtained as an average over the entire ensemble of matrices. We note that we have employed precisely the same data set in performing ensemble and spectral unfolding calculations. Thus, the resulting differences provide a fair measure of differences in the unfolding procedure.

The smoothed spectral density $`\rho (E)`$ is not uniquely defined, and one can think of many methods of equal a priori merit to obtain it. In our case, we have chosen polynomial unfolding. We unfold with a polynomial of fifth order. The order is kept fixed throughout the analysis. Given this unfolding method, an additional parameter must be fixed, namely the length $`L_{\mathrm{fit}}`$ of the interval in which we unfold or, equivalently, the fraction $`L_{\mathrm{fit}}/N`$.

In choosing $`L_{\mathrm{fit}}/N`$, one must be aware of the massive finite sample size effects that can be encountered. In a sample of $`N`$ levels, the number variance vanishes trivially for $`L=N`$. The onset of this effect, however, can already be seen at intervals of length $`LN`$. Numerically studies suggest that interval lengths of $`N/10`$ or less are required for the GOE in order to avoid finite size effects, which cause the number variance to decrease beyond a certain $`L`$. By unfolding in a finite interval, $`L_{\mathrm{fit}}`$, the number of eigenvalues is rather $`L_{\mathrm{fit}}`$ than $`N`$. Unwanted finite size effects can be observed if $`L_{\mathrm{fit}}`$ is chosen too small. Such effects are trivial and do not imply a violation of spectral ergodicity.

Keeping finite size effects in mind, we proceed to investigate the effect of $`L_{\mathrm{fit}}`$ on the spectral statistics. Fig. 8 shows the number variance for fixed $`N=1000`$ and $`d_{\mathrm{eff}}=10`$ as a function of $`L_{\mathrm{fit}}`$. For large intervals, $`L_{\mathrm{fit}}=0.9N`$, the number variance is similar to that obtained from ensemble unfolding. Fluctuations are of Wigner-Dyson form for $`L10`$ and then become stronger. As $`L_{\mathrm{fit}}`$ is reduced from $`0.9N`$ to roughly $`0.3N`$, the onset of the discrepancy between the empirical results and the Wigner-Dyson is displaced to larger values of $`L`$. This is clear from the plot for $`L_{\mathrm{fit}}/N=0.3`$ shown in Fig. 8. Note, however, the appearance of a saturation of the fluctuations below the GOE for large values of $`L`$. This saturation is even more pronounced for $`L_{\mathrm{fit}}/N=0.1`$. The fluctuations in the case $`L_{\mathrm{fit}}/N=0.1`$ for $`L>10`$ are clearly distorted by finite size effects. However, the observation remains that deviations from the GOE are reduced when the length of the fitting interval is reduced moderately. The dramatic differences between the number variance shown here and in the preceding sections are due strictly to differences in the unfolding procedures.

A qualitative explanation of the results of Fig. 8 can be obtained as follows. Fig. 9 shows the difference between the fitted (i.e. smoothed) and integrated spectral densities as a function of the unfolding interval for a single matrix. Standard GOE results are also shown. While fluctuations are obviously present, there is also evidence of strong correlations on all energy scales. Correlations are also clearly present in the GOE results. The structure seen on a macroscopic scale is unwanted. It is in no sense universal and changes markedly from matrix to matrix. The goal of unfolding is thus to remove this non-universal macroscopic structure while preserving the universal microscopic structure of interest. This is evidently a delicate task. If $`L_{\mathrm{fit}}/N`$ is too large, as for $`L_{\mathrm{fit}}/N=0.9`$, a fixed-order polynomial fit is incapable of removing all macroscopic structure. The remaining structure will lead to a number variance larger than that of the GOE. If $`L_{\mathrm{fit}}/N`$ is too small, as for $`L_{\mathrm{fit}}/N=0.3`$, unfolding will begin to eliminate the genuine correlations which we would like to investigate. The resulting number variance will then be smaller than that of the GOE. Evidently, there should be some choice of $`L_{\mathrm{fit}}/N`$ such that the number variance of the GOE is reproduced essentially exactly.

Although the results of Fig. 9 are related to the Thouless energy, they do not provide us with a useful tool for its extraction. They do provide a qualitative explanation of the differences arising from the different unfolding approaches. Ensemble unfolding naturally includes the effects of fluctuations on all scales and will tend to maximize the disagreement with RMT. By contrast, self-unfolding removes long wave length correlations and gives results for the number variance which are in better agreement with random matrix theory. Given its extreme sensitivity to the unfolding procedure, the number variance cannot be offered as evidence for the violation of spectral ergodicity. The safest conclusion is that the choice of unfolding procedure has important consequences and must be made on the basis of physical considerations.

## VII Normal mode analysis

We have seen in the previous sections that ensemble and spectral unfolding do not necessarily lead to the same answers. Moreover, the results of the latter depend on the way in which smoothing of the spectral density is performed. The differences can be understood from the collective motion of large numbers of eigenvalues. The GOE spectrum is remarkably rigid with regard to correlated fluctuations of all wave lengths. While the spectra of sparse random matrices show a similar stiffness to short wave length fluctuations, they are far more susceptible to long wave length fluctuations. The process of self-unfolding involves the elimination of precisely these long wave length modes. Once they have been eliminated, the unfolded spectrum which remains displays a rigidity much closer to that of the GOE. Initial investigations along this line were carried out by French, who associated apparent non-ergodic behavior with the collective motion of eigenvalues. He also found a semi-empirical formula which connects ensemble and spectral averaged number variance . These findings were recently revisited in Ref. .

For our purposes, it is most convenient to consider the set of “normal modes” describing the statistically independent fluctuations of the eigenvalues of a random matrix about their ensemble averaged locations. They provide a complete set of functions suitable for the describing of collective spectral motion . We briefly outline relevant considerations of Ref.. The average positions of the eigenvalues on the unfolded scale is, by construction,

$$\overline{x_i}=i.$$
(19)
The correlation matrix

$$D_{ij}=\overline{x_ix_j}\overline{x_i}\overline{x_j},$$
(20)
provides a measure of the fluctuations of the eigenvalues around their average position. Evidently, the matrix $`D`$ is only easy to define from an ensemble average. Since $`D`$ is a Hermitian matrix, it has $`N`$ real eigenvalues, $`\omega _k`$, obtained from

$$D_{ij}\psi _j^{(k)}=\omega _k\psi _i^{(k)},$$
(21)
with corresponding eigenfunctions $`\psi _i^{(k)}`$. The eigenvalues $`\omega _k`$ measure the mean square amplitude of the corresponding fluctuation and are, of course, quite distinct from the eigenvalues, $`x_i`$, of the random matrix itself. By definition, the $`\psi _i^{(k)}`$ are statistically independent. The actual positions of the eigenvalues of any given matrix in the ensemble can always be written as the sum of their average positions and their fluctuations, $`\delta x_i`$, around them. The fluctuations can be expanded in terms of the eigenfunctions of $`D`$, which gives

$`x_i`$ $`=`$ $`\overline{x_i}+\delta x_i`$ (22)
$`=`$ $`i+{\displaystyle \underset{k=1}{\overset{N}{}}}c_k\psi _i^{(k)},`$ (23)
with some coefficients $`c_k`$. Given the completeness of the $`\psi _i^{(k)}`$, there is a unique correspondence between the $`x_i`$ and the $`c_k`$ for any given matrix. It is clear that $`\overline{c_k}=0`$. Hence, we can express the ensemble number variance Eq. (14) as

$`\overline{\mathrm{\Sigma }}_{L_0}^2(L)`$ $`=`$ $`\overline{(x_{L_0+L/2}x_{L_0L/2})^2}L^2`$ (24)
$``$ $`{\displaystyle \underset{k=1}{\overset{N}{}}}\omega _k\left(\psi _{L_0+L/2}^{(k)}\psi _{L_0L/2}^{(k)}\right)^2,`$ (25)
where we have made use of the fact that

$$\overline{c_kc_k^{}}=\omega _k\delta _{kk^{}}.$$
(26)
We have made a large $`N`$ approximation in obtaining the second form in Eq. (25) and have assumed that the normal modes, $`\psi _i^{(k)}`$ are smooth functions of $`i/N`$. These approximations lead to acceptable errors of size $`1/N`$ in $`\overline{\mathrm{\Sigma }}_{L_0}^2(L)`$.

For the Poisson case with a uniform distribution in the interval $`[0:N+1]`$, Eq. (21) can be solved exactly using the methods introduced in Ref.  <sup>*</sup><sup>*</sup>*This differs slightly from Ref. , where Eq. (21) was solved with the Poisson joint probability density of the eigenvalues as the starting point.. The eigenvalues are

$$\omega _k=\frac{N+1}{4(N+2)\mathrm{sin}^2\left(\pi k/2(N+1)\right)}$$
(27)
and the corresponding normal modes are

$$\psi _i^{(k)}=\sqrt{\frac{2}{N+1}}\mathrm{sin}\left(\frac{\pi ki}{N+1}\right).$$
(28)
For the GOE, it is convenient to adopt a slightly different approach. Instead of constructing the matrix $`D`$ of Eq. (21), it is easier to start from the matrix

$$C_{ij}=\frac{^2}{x_ix_j}\mathrm{log}P_N(x_1,\mathrm{},x_N),$$
(29)
where $`P_N(x_1,\mathrm{},x_N)`$ is the joint probability distribution of the eigenvalues. The $`C`$-matrix describes the fluctuations of the eigenvalues around their average positions. The corresponding eigenvalue equation of $`C`$ can be solved exactly. The matrices $`C`$ and $`D`$ are closely related. To the extent that the small amplitude approximation is valid, the eigenfunctions of the matrices $`C`$ and $`D`$ are identical; the corresponding eigenvalues are simply reciprocals. Thus, although $`C`$ and $`D`$ are qualitatively similar, they do not necessarily lead to identical answers. The eigenvalues of $`C`$ for the GOE are thus given by

$$\omega _k=\frac{N}{2\sqrt{2}\pi k},$$
(30)
which is slightly different from the corresponding result in Ref.  due to our use of a different mean level spacing at $`E=0`$. Taking constant factors into account, the normal modes of the GOE eigenvalues are given in the large $`N`$ limit by

$$\psi _i^{(k)}=\sqrt{\frac{2}{N}}U_{k1}\left(\frac{\pi (iN/2)}{N}\right),$$
(31)
where $`U_n(x)`$ are the Chebyshev polynomials of the second kind.

With these explicit expressions, the qualitative interpretation of the normal modes becomes obvious. The form of the eigenfunctions suggests that we can roughly associate a wave length, $`1/k`$, with each normal mode. The normal modes then appear as compressional waves in the spectra. In each case, the mode of longest wave length (i.e. with $`k=1`$) involves all eigenvalues moving in the same direction. The eigenvalues, $`\omega _k`$, measure the mean square amplitude of the fluctuations of the various normal modes. It is clear from Eq. (25) that the largest contributions to the number variance will come from modes with the largest mean square amplitude. In the case of the GOE, this mean square amplitude is strictly proportional to $`1/k`$ and provides striking confirmation of the “rigidity” associated with the spectra of the Gaussian ensembles. It is precisely this feature of the spectrum of normal modes which gives rise to the logarithmic asymptotic behavior of the number variance. By contrast, the long wave length normal modes of the Poisson distribution indicate $`\omega _kN^2/\pi ^2k^2`$. Eq. (25) now makes it clear that this softness of long wave length modes is directly responsible for the linear asymptotic behavior of the number variance. In short, information regarding the fluctuations of eigenvalues about their average positions as described by the normal modes provides a compact source of information regarding longer range spectral fluctuations.

We now turn to the numerical solution of Eq. (21) with our present data and a comparison with the theoretical predictions. Fig. 10 shows the resulting dispersion relations for both the GOE and Poisson cases. In each case, the lowest 100 modes are shown for several matrix sizes. The agreement between theory and numerical data is quite impressive. In the case of the GOE, deviations from the analytical predictions for small matrix dimensions are due to the limitations of the Gaussian approximation.

These results lead immediately to the important observation that the eigenfunctions of $`D`$ for arbitrary $`d_{\mathrm{eff}}`$ are reasonably well described by the sine waves of Eq. (28) independent of whether the spectral statistics are pure GOE, pure Poisson, or a scale-dependent mixture of the two. As noted, the eigenvalues, $`\omega _k`$, depend sensitively on the statistics chosen (i.e. on $`d_{\mathrm{eff}}`$). It is therefore reasonable to approximate the number variance at $`L_0=0`$, according to Eq. (25), as

$$\overline{\mathrm{\Sigma }}^2(L)\frac{4}{N}\underset{k=1}{\overset{N}{}}\omega _k\mathrm{sin}^2\left(\frac{\pi Lk}{2N}\right).$$
(33)
In this form, we see that the number variance is largely determined by the dispersion relation of the eigenvalues $`\omega _k`$. In particular, the large $`L`$ behavior is dominantly set by the softest normal modes.

The normal mode spectra and Eq. (33) can now provide us with a simple explanation of the differences between the number variance obtained from ensemble and spectral unfolding in the case of matrices with sparsity $`0<d_{\mathrm{eff}}N`$. Let us first consider the case of ensemble unfolding. First, consider $`L`$ to be of order $`1`$. The contribution of long wave length modes, for which $`\omega _k1/k^2`$, is suppressed by the factor $`k^2`$ coming from the sine term. The contribution from short wave length modes (i.e. $`k`$ of order $`N`$), for which $`\omega _k1/k`$, persists and builds the logarithmic behavior of the number variance familiar from the GOE. As $`L`$ increases, the long wave length modes are no longer suppressed. Their $`1/k^2`$ contributions now build up the linear divergence of the number variance familiar from the Poisson distribution. The resulting qualitative behavior is precisely that shown in Fig. 4.

In Fig. 14, we give some examples of the number variance obtained from evaluating the sum Eq. (33). We consider the GOE, Poisson, and two specific toy dispersion relations motivated by the results of our sparse matrix model. The first of these is intended to mimic the results of ensemble unfolding. We construct a dispersion relation which interpolates between GOE and Poisson forms, i.e.

$$\omega _k=\{\begin{array}{cc}\omega _k^{(\mathrm{Poisson})}/k_\mathrm{c}\hfill & kk_\mathrm{c}\\ \omega _k^{(\mathrm{GOE})}\hfill & k>k_\mathrm{c}\end{array},$$
(34)
where $`\omega _k^{(\mathrm{Poisson})}`$ and $`\omega _k^{(\mathrm{GOE})}`$ are given by Eqs. (27) and (30), respectively. The choice $`k_\mathrm{c}=\sqrt{N}`$ leads to a nearly continuous dispersion relation and a surprisingly faithful reproduction of the dispersion relation shown in Fig. 12 for the case $`d_{\mathrm{eff}}=10`$. The resulting number variance is in striking agreement with the ensemble average results shown in Fig. 4.

The second toy dispersion relation is intended to demonstrate the striking sensitivity of the number variance to the treatment of the very softest normal modes through the introduction of a simple cutoff

$$\omega _k=\{\begin{array}{cc}0\hfill & kk_0\\ \omega _k^{(\mathrm{GOE})}\hfill & k>k_0\end{array}.$$
(35)
It is evident from Eq. (33) that every normal mode makes a positive contribution to $`\overline{\mathrm{\Sigma }}^2(L)`$. The elimination of soft modes resulting from the use of Eq. (35) will necessarily lead to a number variance everywhere smaller than the GOE result. Even the modest choice of $`k_0=10`$ for $`N=1000`$ shown in Fig. (14) is sufficient to cause the number variance to saturate and decrease for sufficiently large $`L`$.

These simple examples lead us towards a better understanding of the results of self-unfolding. As we have seen, it is necessary to eliminate non-universal, long wave length fluctuations in the spectrum (i.e. unfold the spectrum) if we are to be able to make sensible comparisons of spectral fluctuations in macroscopically distinct regions of the spectrum. Let us now consider unfolding from the point of view of Eq. (23). The eigenvalues from any given realization of a random matrix uniquely determine the expansion coefficients, $`c_k`$, in the complete set of normal modes. Self-unfolding then corresponds to imposing a smooth cutoff on these coefficients to eliminate the contribution of normal modes with small $`k`$. In other words, one replaces $`c_k`$ by $`G(k)c_k`$ where $`G(k)0`$ as $`k0`$. Obviously, the scale on which $`G(k)`$ vanishes must be set by physical arguments. If desired, of course, one can invert Eq. (23) and obtain the resulting unfolded eigenvalues. Qualitatively, one can proceed directly to the evaluation of the number variance through the approximate Eq. (33). Using Eq. (26), it is clear that the primary effect of the unfolding just described is to replace the terms $`\omega _k`$ by $`\omega _kG^2(k)`$. This has the effect of greatly reducing the contribution of long wave length modes to the number variance. However, these are precisely the modes which are responsible for the large $`L`$ form of the number variance.

There are now several possibilities. If the normal mode spectrum is scale invariant, the suppression of long wave length modes implied by unfolding will have little consequence. The number variance will have the same qualitative behavior independent of the unfolding method adopted. This is the case for both the Poisson distribution and the GOE. The situation is quite different for sparse random matrices where the normal mode spectrum displays a transition from Poisson to Gaussian form on a scale $`k_c\sqrt{N}`$. The number variance obtained from ensemble unfolding will show deviations from random matrix theory. If the scale of $`G(k)`$ implicit in self-unfolding is sufficiently large to eliminate the contributions of soft normal modes, we will obtain the logarithmic behavior of the Gaussian ensembles. While such disagreement has been viewed as a breaking of spectral ergodicity, this is not the case. Rather, ensemble and spectral unfolding explore the same normal mode spectra at different scales. If different questions are asked in the two cases, we should not be surprised to get different answers. Although it is conventional to view self-unfolding as the more ambiguous procedure, it is also possible to restore agreement between the two methods by modifying the ensemble unfolding procedure. Specifically, this could be accomplished by increasing the bin size used in Eq. (12) so that each bin included roughly $`\sqrt{N}`$ levels.

The effect of polynomial self-unfolding on the normal modes can be seen in Fig. 15. Here, we construct the matrix $`D`$ according to Eq. (20) using the self-unfolded eigenvalues. The normal modes in Fig. 15 are obtained from the same data used in constructing Fig. 12. The effects are quite dramatic but entirely predictable. The purpose of unfolding is to eliminate or greatly reduce the mean square amplitude of long wave length normal modes (i.e. reduce $`\omega _k`$ for small $`k`$.) This figure reveals that the softest modes are lowered by one order of magnitude for an unfolding interval of length $`L_{\mathrm{fit}}=0.9N`$. For $`L_{\mathrm{fit}}`$ of $`0.4N`$ to $`0.5N`$, the relation dispersion is remarkably similar to that of the GOE, which is consistent with our results for $`\mathrm{\Sigma }^2`$. For even smaller intervals, the dispersion relation falls below the GOE and gives rise to the saturation effects seen in Fig. 8. It is clear that the effective cutoff, $`G(k)`$ can be determined immediately as the ratio of the results of Fig. 15 to those of Fig. 12.

## VIII Summary and conclusions

We have presented a detailed analysis of the spectral properties of sparse random matrix ensembles with particular emphasis on the spectral ergodicity hypothesis and the effects of sparsity on spectral statistics. In addition, we have presented the first numerical investigation of a normal mode spectrum in random matrix ensembles. The question of spectral ergodicity is of direct relevance for a variety of physical systems. The existence of apparent counter examples, such as lattice QCD and many-body systems with two-body interactions, suggests that this hypothesis is not always fulfilled. The correct averaging procedure seems clear in these two cases. In lattice QCD, the ensemble average is the relevant procedure by construction. Spectral averaging seems physically meaningless, and differences between the approaches have not been regarded as surprising. On the other hand, the spectral averaging is evidently more appropriate for the TBRE, since we are interested in an “average nuclear spectrum” and not a spectrum averaged over many distinct nuclei. In each of these cases, the physically motivated averaging procedure agrees with data. In the case of disordered systems, the ergodicity hypothesis is a crucial ingredient necessary for the comparison of theory and experiment. Due to the similarities of all three systems, e.g. with respect to spectral statistics, it is important to understand if disordered systems can show non-ergodic properties or if, on the other the hand, TBRE and lattice QCD theory are in reality ergodic. Our analysis indicates the latter and suggests that apparent differences between spectral and ensemble averaged results have been misinterpreted.

Using ensemble unfolding, we presented a detailed analysis of the parameter dependence of our model. For a certain parameter range, spectral correlations are described by RMT; outside this range, fluctuations are stronger. The dependence of this critical energy scale on the model parameters can be understood from viewing the model as a multi-dimensional disordered system on a random lattice with fixed disorder strength. Results obtained from a self-unfolding of the spectra do not agree with the results of ensemble unfolding. There is also a strong dependence on the parameters of the self-unfolding procedure, i.e. the length of the unfolding interval. Changes of this length can alter the critical energy obtained with ensemble unfolding. Indeed, this length can even be chosen to reinstate virtually perfect agreement with the spectral correlations of random matrix theory.

The observed differences between the results of ensemble and spectral unfolding have been shown to be a result of the collective motion of large numbers of eigenvalues. Such differences are intrinsic to systems lacking scale invariance, for which correlations change from RMT to Poisson on some scale. Ensemble unfolding includes the effects of fluctuations on all scales. Self-unfolding necessarily eliminates the effects of “soft”, long wave length fluctuations, whose contributions to the asymptotic number variance would otherwise be dominant. We have shown that the results of self-unfolding are likely to be closer to those of RMT and are likely to differ from the results of ensemble unfolding. Spectral ergodicity is not broken. Rather, the two approaches probe different scales of the spectrum.

The normal modes of the random matrix proved to be a convenient tool for the investigation of spectral correlations. They describe the fluctuations of the eigenvalues around their average positions and can be understood approximately as compressional waves. Our analysis shows that the normal modes are well described as plane waves, independent of the parameters of our model. Their dispersion relation, however, is sensitive to the choice of parameters. Spectral correlations are therefore largely determined by dispersion relation, as indicated by the explicit expression for the number variance given above. The presence of a clear scale dependence in the dispersion relation for sparse random matrices was sufficient to provide a qualitative explanation of the apparent violation of spectral ergodicity in this problem. This example serves a striking reminder that the unfolding procedure, usually regarded as technically difficult and uninteresting, can have important physical content. An appropriate unfolding procedure should reflect the spectral scale which is relevant for the physical properties in question. This scale is not always apparent given the usual ad hoc treatment of unfolding. The normal modes permit a more systematic formulation of this problem.

Given the similarities between disordered systems and the sparse random matrices considered here and the matrices appropriate for disordered systems, it would be of interest to perform a normal mode analysis appropriate for this case as well. A scale invariant dispersion relation could provide strong confirmation of spectral ergodicity in disordered systems; the presence of a scale could indicate interesting new lines of experimental investigation.

## Acknowledgments

We are grateful to acknowledge discussions with T. Seligman and T. Papenbrock.


format: text

text:
# References

Reply to Comment on Dirac spectral sum rules for QCD in three dimensions

U. Magnea

INFN, Via P. Giuria 1, I-10125 Torino, Italy
blom@to.infn.it

I reply to the comment by Dr S. Nishigaki (hep-th/0007042) to my papers Phys. Rev. D61 (2000) 056005 and Phys. Rev. D62 (2000) 016005.

Recently, Nishigaki  pointed out that the first sum rule for the eigenvalues of the Dirac operator given in my papers  was incorrect. This appears to be true, and I thank Dr Nishigaki for drawing attention to this fact. That the sum rules in  do not agree with the results by Hilmoine and Niclasen  had been communicated to me by these authors already before the appearance of their paper, but I was at the time unable to detect where my calculation failed.

However, I do not agree with the explantion by Dr Nishigaki of why the sum rules were not reproduced correctly. According to him, the authors of  “were aware that they needed to extend the saddle point manifold… in order to obtain a correct result”. Such an extension of the integral over the coset obviously can not have any effect on the result, since this procedure only involves integrating over dummy variables on which the integrand does not depend (namely the degrees of freedom of the stability subgroup), and amounts only to multiplying the integral by a numerical factor. This extension was done in  only for convenience. Also, to my understanding, Dr Nishigaki’s references to the papers by Kamenev and Mèzard are not relevant in this context.

Instead, the reason is quite simple. In , formula (6.47), and in , formula (8.60), in order to obtain the right sum rule, the sum should go from $`1`$ to $`M_{as}`$ and from $`1`$ to $`M_s`$ respectively, where $`M_{as}`$ is the number of antisymmetric traceless generators that generate antisymmetric unimodular matrices, and $`M_s`$ is the number of symmetric traceless generators that generate symmetric unimodular matrices. Taking into account the sizes of these matrices in the two cases we get

$$M_{as}=\frac{4N_f(4N_f1)}{2}1,M_s=\frac{2N_f(2N_f+1)}{2}1$$
(1)
We then obtain for $`\beta =1`$

$$\underset{\lambda _k>0}{}\frac{1}{(N\mathrm{\Sigma }\lambda _k)^2}=\frac{2N_f}{2(2N_f1)(4N_f+1)}$$
(2)
and for $`\beta =4`$

$$\underset{\lambda _k>0}{}\frac{1}{(N\mathrm{\Sigma }\lambda _k)^2}=\frac{2N_f}{2(2N_f1)(N_f+1)}$$
(3)
This is in total agreement with the sum rules obtained by Hilmoine and Niclasen and by Nishigaki. I remind the reader that my definition of $`N_f`$ differs by a factor of two from the other authors’ definition, and that I sum only over the positive eigenvalues. In the sum rule for $`\beta =1`$, in addition to the above mentioned correction, the sum rule was off by a factor of 4 due to various mistakes of factors of 2 for which I apologize. I will update the electronic versions of  shortly.

The $`3d`$ sum rules can be summarized as follows:

$$\underset{\lambda _k>0}{}\frac{1}{(N\mathrm{\Sigma }\lambda _k)^2}=\frac{2N_f}{2(2N_f1)(\frac{4N_f}{\beta }+1)}$$
(4)
where the number of flavors is always even and denoted by $`2N_f`$, and $`\beta `$ is the Dyson index of the corresponding matrix model. The sum rule for $`\beta =2`$ was given by Verbaarschot and Zahed in .

Acknowledgment

I thank Jac Verbaarschot for correspondence. I am indebted to him for assisting me in finding the mistake in my sum rules .


text:
# Exact particle and kinetic energy densities for one-dimensional confined gases of non-interacting fermions

## I General formulation

The one-body Dirac density matrix for a system of $`N`$ non-interacting fermions at zero temperature can be expanded on the single-particle wavefunctions $`\psi _i(x)=x|\psi _i`$ as $`\rho (x_1,x_2)=_{i=1}^N\psi _i^{}(x_1)\psi _i(x_2)`$. By using the representation of the translation operator this becomes

$$\rho (x_1,x_2)=\underset{i=1}{\overset{N}{}}\psi _i^{}(x_1)e^{i\widehat{p}(x_1x_2)}\psi _i(x_1),$$
(1)
showing how distant points are correlated through the momentum operator $`\widehat{p}`$. Expansion in powers of the relative coordinate $`r=x_1x_2`$ yields physical observables such as the particle density profile $`n(x)`$,

$$n(x)=\rho (x+r/2,xr/2)|_{r=0}=\underset{i=1}{\overset{N}{}}\psi _i|\delta (xx_i)|\psi _i$$
(2)
and the kinetic pressure $`P(x)`$,

$$P(x)=\frac{\mathrm{}^2}{m}\frac{^2}{r^2}\rho (x+r/2,xr/2)|_{r=0}=\frac{1}{2m}\underset{i=1}{\overset{N}{}}\psi _i|p_i^2\delta (xx_i)+\delta (xx_i)p_i^2|\psi _i.$$
(3)
This is twice the kinetic energy density.

The main idea of this Letter is to rewrite Eqs. (2) and (3) as the imaginary part of the ground-state average of suitable operators related to the Green’s function in coordinate space $`G(x)=(x\widehat{x}+i\epsilon )^1`$. We have

$$n(x)=\frac{1}{\pi }\underset{\epsilon 0^+}{lim}\mathrm{Im}\underset{i=1}{\overset{N}{}}\psi _i|G(x)|\psi _i$$
(4)
and

$$P(x)=\frac{1}{\pi }\underset{\epsilon 0^+}{lim}\mathrm{Im}\underset{i=1}{\overset{N}{}}\psi _i|\widehat{\frac{p^2}{m}}G(x)|\psi _i.$$
(5)
$`G(x)`$ can then be treated by methods analogous to those used for treating Green’s functions in the energy domain.

The equivalence between expressions (2) and (4) is easily proved in the coordinate representation, where the density profile in Eq. (4) reads

$$n(x)=\frac{1}{\pi }\underset{\epsilon 0^+}{lim}\mathrm{Im}\underset{i=1}{\overset{N}{}}dx_i|\psi _i(x_i)|^2\frac{1}{xx_i+i\epsilon },$$
(6)
yielding Eq. (2) when one takes the limit $`\epsilon 0^+`$. The equivalence between expressions (3) and (5) for $`P(x)`$ is similarly proved.

Evidently, this method can be applied to all 1D systems which may be described by single-particle orbitals: one only needs to know the representation of the position and momentum operators on such a basis. Hence, interactions could also be included in evaluating the particle density through the use of Kohn-Sham single-particle orbitals. Models of displacement fields (such as those induced by impurities) may also be studied directly without previous evaluation of orbitals.

## II Non-interacting Fermi gas in harmonic trap

As already noted, a 1D Fermi model is relevant to the spin-polarized fermionic vapors in magnetic confinement , where it is possible to realize experimentally a 1D configuration by making use of very anisotropic axially symmetric traps. At low temperature only the transverse ground state of the trap is populated and the vapor can be described by an effective 1D harmonic Hamiltonian.

An analytic expression for the particle density of this system has been given by Husimi  in terms of the wavefunction of the $`N`$-th fermion in the trap. However, a calculation of the density profile and the kinetic pressure by his approach is limited to small values of $`N`$ and has been reported only for $`N=`$ 1 and 2 . Our method allows us to efficiently evaluate these ground state properties even for quite large numbers of particles. Of course, the simplicity of the representation of the position and momentum operators in this system makes it a favorable example.

### A Particle density profile

We start with the calculation of the particle density. For a linear harmonic oscillator the position operator is represented as $`\widehat{x}=(a+a^{})/\sqrt{2}`$ on the basis $`\{|\psi _i\}`$ of the eigenvectors of the Hamiltonian. As usual, the creation and destruction operators satisfy the relations $`a|\psi _i=\sqrt{i1}|\psi _{i1}`$ and $`a^{}|\psi _i=\sqrt{i}|\psi _{i+1}`$. A straightforward procedure for evaluating the profile (4) is to invert the matrix $`(x\widehat{x}+i\epsilon )`$ and to calculate its trace on the submatrix of its first $`N\times N`$ block (TrN). We employ the relation $`\mathrm{Tr}_NQ=\left[\mathrm{ln}det(Q^1+\lambda 𝕀_N)\right]/\lambda |_{\lambda =0}`$, where $`𝕀_N`$ is a diagonal semi-infinite matrix with its first $`N`$ eigenvalues equal to 1 and null elsewhere , to obtain

$$n(x)=\frac{1}{\pi }\underset{\epsilon 0^+}{lim}\mathrm{Im}\frac{}{\lambda }\left[\mathrm{ln}det(x\widehat{x}+i\epsilon +\lambda 𝕀_N)\right]_{\lambda =0}.$$
(7)
The calculation of the determinant in Eq. (7) is conveniently performed by the recursive algorithm developed in . Renormalization of the tridiagonal operator $`\widehat{R}=\widehat{x}\lambda 𝕀_N`$ allows us to write

$$det(x\widehat{R}+i\epsilon )=\underset{k=1}{\overset{\mathrm{}}{}}(x\stackrel{~}{a}_k+i\epsilon )$$
(8)
with $`\stackrel{~}{a}_1=\lambda `$, $`\stackrel{~}{a}_{k+1}=\lambda +\frac{1}{2}k(x\stackrel{~}{a}_k+i\epsilon )^1`$ for $`1<k<N`$ and $`\stackrel{~}{a}_{k+1}=\frac{1}{2}k(x\stackrel{~}{a}_k+i\epsilon )^1`$ for $`kN`$.

The scheme given in Eqs. (7) and (8) is easily implemented numerically. In practice we have performed the calculation of the determinant (8) up to the product of its first $`M`$ terms, which corresponds to inverting an $`M\times M`$ matrix. We have checked the convergence of this approximation by increasing the dimension $`M`$ and correspondingly decreasing the value of $`\epsilon `$.

In Fig. 1 we report the density profile $`n(x)`$ for $`N=5`$, 10 and 20 fermions, with $`M=10^5`$ and $`\epsilon =0.01`$. The exact profiles are also compared with those given by the Thomas-Fermi approximation (LDA),

$$n_{\mathrm{LDA}}(x)=\frac{1}{\pi }(2Nx^2)^{1/2}$$
(9)
(in units such that $`\mathrm{}=1`$, $`m=1`$ and $`\omega =1`$). The exact profiles contain $`N`$ oscillations, which become smaller in relative amplitude as $`N`$ increases. Without any special numerical efforts we have evaluated the exact density profile up to $`N=1000`$: this would otherwise require the calculation of Hermite polynomials up to the 1000th degree. For large $`N`$ the oscillations are so small in relative amplitude that their smoothing in the LDA profile becomes reasonably accurate, except for the region of the tails.

### B Kinetic pressure profile

The particle density profile that we have evaluated above is the analogue of the density of single-particle states in the energy domain. In the space domain other single-particle quantities acquire physical interest, as is the case for the kinetic pressure in Eq. (3). We show how also for this function one can profitably resort to a renormalization technique.

Taking $`Q=\widehat{p}^2G(x)`$, Eq. (5) can be written as

$$P(x)=\frac{1}{\pi }\underset{\epsilon 0^+}{lim}\mathrm{Im}\frac{}{\lambda }\left[\mathrm{ln}det(x\widehat{x}+i\epsilon +\lambda 𝕀_N\widehat{p}^2)\right]_{\lambda =0}.$$
(10)
The matrix $`(x\widehat{x}+i\epsilon +\lambda 𝕀_N\widehat{p}^2)`$ appearing in Eq. (10) is pentadiagonal on the first $`N`$ rows, owing to the form of the operators $`\widehat{x}`$ and $`\widehat{p}=i(a^{}a)/\sqrt{2}`$ in the basis of the energy eigenstates $`\{|\psi _i\}`$. The calculation of the determinant of such a matrix can again be performed by the recursive algorithm given in . Renormalization of the operator $`\widehat{K}=\widehat{x}\lambda 𝕀_N\widehat{p}^2`$ is made on blocks of dimension 2 for the pentadiagonal part and on blocks of dimension 1 for the tridiagonal part. This allows us to write

$$det(x\widehat{K}+i\epsilon )=\{\begin{array}{c}_{j=1}^{(N+2)/2}det(x\stackrel{~}{A}_j+i\epsilon )_{k=N+3}^{\mathrm{}}(x\stackrel{~}{a}_k+i\epsilon )(\mathrm{even}N)\hfill \\ _{j=1}^{(N+3)/2}det(x\stackrel{~}{A}_j+i\epsilon )_{k=N+4}^{\mathrm{}}(x\stackrel{~}{a}_k+i\epsilon )(\mathrm{odd}N).\hfill \end{array}$$
(11)
The renormalized $`2\times 2`$ blocks $`\stackrel{~}{A_j}`$ satisfy the recursion relation

$$\stackrel{~}{A_j}=A_j+B_{j,j1}(x\stackrel{~}{A}_{j1}+i\epsilon )^1B_{j1,j}$$
(12)
for $`j>1`$ and $`\stackrel{~}{A}_1=A_1`$. The matrices $`A_j`$, $`B_{j,j+1}`$ and $`B_{j+1,j}`$ are submatrices of the operator $`\widehat{K}`$, which are defined as follows:

$$A_j=\left(\begin{array}{cc}\lambda \left(2j3/2\right)\theta _{N2j+1}& \sqrt{(2j1)/2}\\ \sqrt{(2j1)/2}& \lambda \left(2j1/2\right)\theta _{N2j}\end{array}\right),$$
(13)
$$B_{j,j+1}=\left(\begin{array}{cc}\lambda \sqrt{2j(2j1)}/2\theta _{N2j+1}& 0\\ \sqrt{j}& \lambda \sqrt{2j(2j+1)}/2\theta _{N2j}\end{array}\right)$$
(14)
and

$$B_{j+1,j}=\left(\begin{array}{cc}\lambda \sqrt{2j(2j1)}/2\theta _{N2j1}& \sqrt{j}\\ 0& \lambda \sqrt{2j(2j+1)}/2\theta _{N2j2}\end{array}\right),$$
(15)
with $`\theta _k=1`$ for $`k0`$ and $`\theta _k=0`$ otherwise. The recursion relation for the elements $`\stackrel{~}{a}_k`$ is again $`\stackrel{~}{a}_{k+1}=\frac{1}{2}k/(x\stackrel{~}{a}_k+i\epsilon )`$, the first elements being $`\stackrel{~}{a}_{k+1}=\frac{1}{2}k\{x+i\epsilon [\stackrel{~}{A}_{k/2}]_{22}[\stackrel{~}{A}_{k/2}]_{21}[\stackrel{~}{A}_{k/2}]_{12}/(x+i\epsilon [\stackrel{~}{A}_{k/2}]_{11})\}^1`$ with $`k=N+2`$ for even $`N`$ and $`k=N+3`$ for odd $`N`$. We have studied the convergence of the determinant in Eq. (11) as for the case of the particle density profile.

In Fig. 2 the kinetic pressure $`P(x)`$ is plotted for $`N=5`$, 10 and 20 with $`M=10^5`$ and $`\epsilon =0.01`$, together with the profiles $`P_{\mathrm{LDA}}(x)`$ evaluated in the Thomas-Fermi approximation,

$$P_{\mathrm{LDA}}=\frac{1}{3\pi }(2Nx^2)^{3/2}.$$
(16)
The exact kinetic pressure shows $`N`$ oscillations and has the peculiarity of being negative in the tails. This microscopic quantum effect, which is missing in the local density description, reflects the fact that in the low density region the kinetic energy decreases with increasing density. We have checked that our results agree with those reported in  for $`N=1`$ and 2, and carried out the calculation of $`P(x)`$ up to $`N=1000`$. The kinetic pressure profile for N=1000, as shown in Fig. 3, looks almost indistinguishable from the LDA prediction but stills presents a region of negative kinetic pressure in the tails (inset in Fig. 3).

In conclusion, in this Letter we have given a general formula for the exact particle density and kinetic pressure profiles of a 1D many-fermions system in terms of a Green’s operator in coordinate space. We have made use of the decimation/renormalization procedure and of other recursive techniques, originally developed to evaluate the spectral properties of quasi-1D systems in solid state physics, to efficiently calculate the exact density profiles of a harmonically confined non-interacting Fermi gas. Within the same general scheme the particle density could also be evaluated by employing a suitable Kirkman-Pendry relation , as will be reported elsewhere. We have verified that for large number of atoms ($`N=1000`$) the local density approximation reproduces reasonably well the exact profiles except for the region of the tails, where the exact kinetic pressure is negative.

We believe that the present method opens the way for a novel approach to the equilibrium properties of spatially inhomogeneous 1D systems. The expressions here derived can be extended to finite temperature and to calculate partial density profiles for subgroups of atoms. The kinetic energy density functional can be studied through the calculation of the function $`P[x(n)]/2`$ where $`x(n)`$ is obtained by local inversion of the exact profile $`n(x)`$. Pressure fluctuations will become accessible to study through the evaluation of higher moments of the one-body density matrix. The density profiles of the harmonically trapped Fermi gas in 1D, showing a prominent shell structure as displayed in our calculations, could become observable in experiments on alkali vapours.

###### Acknowledgements.

It is a pleasure to thank Professor N. H. March for useful discussions on Fermi vapors and for drawing our attention to the work published in Refs. . One of us (P. V.) would like to thank Professor G. Grosso and Dr R. Farchioni for fruitful discussions.