Vibrational spectrum of Granular packings with random matrices

Abstract The vibrational spectrum of granular packings can be used as a signature of the jamming transition, with the density of states at zero frequency becoming nonzero at the transition. It has been proposed previously that the vibrational spectrum of granular packings can be approximately obtained from random matrix theory. Here, we show that the autocorrelation function of the density of states shows good agreement between dynamical numerical simulations of frictionless bead packs near the jamming point and the analytic predictions of the Laguerre orthogonal ensemble of random matrices; there is clear disagreement with the Gaussian orthogonal ensemble, establishing that the Laguerre ensemble correctly reproduces the universal statistical properties of jammed granular matter and excluding the Gaussian orthogonal ensemble. We also present a random lattice model which is a physically motivated variant of the random matrix ensemble. Numerical calculations reveal that this model reproduces the known features of the vibrational density of states of frictionless granular matter, while also retaining the correlation structure seen in the Laguerre random matrix theory. Graphic abstract


INTRODUCTION
Granular materials are a class of systems which are out of equilibrium and not easy to understand within the framework of standard statistical mechanics.For static assemblies, the distribution of forces [1] and the continuum limit [2] are difficult to obtain.This is because interparticle contacts are very stiff: a slight compression of two particles that are contact, by an amount that is much less than the interparticle separation, gives rise to large forces.Added complications are caused by the fact that, for noncohesive granular matter, two particles in contact repel each other when they are compressed, but do not attract each other when they are moved away from each other and the contact is broken; that the repulsive force between particles is not a linear function of their compression when the compression is small; [3] and that there are frictional forces between particles, [4] resulting in history dependent forces.The dynamic properties of granular matter are difficult to understand because interparticle collisions are strongly inelastic.If a high density of particles builds up in a region because of random fluctuations, the collision rate and therefore the rate of energy loss increases in the region.This can trap particles in the region, causing the density fluctuations to grow.[5] Experimentally, one observes distinctive phenomena such as force chains and stability against mechanical collapse in very sparse static packings, non-Maxwellian velocity distributions [6][7][8] and inelastic collapse in dilute granular gases, [9,10] and shear thinning and shear thickening in the intermediate regime.[11] The tendency of flowing granular matter to get 'jammed' and stop flowing at low densities is a practical problem that limits the flow rate in the industrial use of granular materials.[12] Remarkably,the transition from a flowing to a jammed state in granular matter, structural glasses, and foams and colloids, can be studied with a unified approach.[13] When the transition occurs at zero temperature and zero shear stress as the density is varied, the transition point is called "Point J", [14] and is characterized by diverging length scales [15,16] suggestive of a second order phase transition.At the same time, other properties of the system change discontinuously at Point J, [14] as one would expect at a first order phase transition.
The density of states for vibrational modes in a granular system is one of the properties that has a signature of the transition at Point J.A jammed granular system has mechanical rigidity.Even though the force between two particles is a nonlinear function of the compression between them, the small deviations from the jammed state (which already has non-zero compression) can be analyzed using a linear model, resulting in normal modes.Extensive numerical simulations [14] on systems at zero temperature and zero shear stress show that the density of states D(ω) as a function of ω approaches zero linearly as ω → 0 if the particle density is greater than the critical density.As the particle density is reduced, the slope of D(ω) at the origin becomes steeper, until at Point J, D(ω → 0) = 0.
In the linearized analyis of vibrational modes, the system can be treated as a network of random springs, with the number of springs decreasing as Point J is approached.It is natural to analyze the problem using ran-arXiv:2006.16497v1[cond-mat.dis-nn]30 Jun 2020 dom matrices, and see how the resultant density of states evolves near the transition.This has been done, [17] and yields a broad peak in the density of states that reaches ω = 0 as the transition is approached.However, the model also predicts a gap in the density of states near ω = 0 above the transition, which does not match the numerical results.There are a few other qualitative discrepancies.
Although it is encouraging that some features of the random matrix density of states matches D(ω) from numerical simulations, it is well known that the overall density of states predicted by random matrix theory often differs from what is observed in systems to which the theory applies, because of non-universal effects [18].Instead, the correlations in the density of states and the distribution of level spacings is a more reliable indicator of the validity of the random matrix approach [18].
In this paper we therefore turn to the correlations in the density of states predicted by random matrix theory.We argue that the Laguerre orthogonal ensemble rather than the Gaussian orthogonal ensemble (GOE) is the appropriate random matrix model for granular bead packs.The correlation function for the Laguerre ensemble differs from that for the GOE near the low-frequency edge of the allowed range of ω [19].By comparing the correlations in the numerically computed vibrational spectrum of granular bead packs near the jamming transition to the predictions of the Laguerre ensemble and the GOE we are able to demonstrate good agreement with the former and to exclude the latter.
The distribution of consecutive level spacings is also a universal feature of the spectrum that should be described by random matrix theory.We find that the level spacing distribution predicted by the Laguerre ensemble is very close to the GOE result both near the zero frequency edge and at high frequency.The spectra calculated for granular bead packs in dynamical simulations are found to agree with this distribution.This finding further validates the random matrix approach but it does not help discriminate between the Laguerre ensemble and the GOE.The agreement of the level spacing distribution with the GOE result has been observed earlier, [20,21] but without reference to the Laguerre ensemble.
We also construct a random lattice model, which is a physically motivated variant of the random matrix ensemble.Although it is not possible to calculate the properties of this model analytically, numerical results reveal that all the qualitative features of D(ω) are reproduced.At the same time, the correlation functions and the level spacing distribution seen in the idealized random matrix theory are not significantly changed.
The rest of the paper is organized as follows.In section II we summarize the random matrix approach to the problem and the results for the density of states.In Section III we compare the autocorrelation function for the density of states and the level spacing distribution for the Laguerre ensemble and the GOE to the vibrational spectra for granular packs near point J, obtained from dynamical numerical simulations.In Section IV we introduce the random lattice model and show that it reproduces both universal and non-universal features of vibrational spectrum of granular packs.In Appendix A we present an analysis of the level spacing distribution for the Laguerre ensemble and in Appendix B some technical details regarding the autocorrelation.

II. LAGUERRE ENSEMBLES
We follow the approach of Ref. [17] here.Within linear response, if the particles in the granular assembly are displaced slightly from their resting positions, their accelerations are of the form ẍ = −Kx, where x is a N component column vector and K is a N × N matrix.The crucial observation [22,23] is that the connection between accelerations and displacements is a two-step process.Within linear response, each contact between a pair of particles can be represented as a spring that has been precompressed by some amount.Thus one has a network of springs, with various spring constants.When a particle is displaced, it stretches (or compresses) each spring that it is connected to, by an amount that is equal to the component of its displacement along that spring.The spring exerts a restoring force that is proportional to this stretching; the spring constant can be different for each spring.Thus we have where Ãij = cos θ ij if the i'th spring is connected to the j'th particle, with θ ij being the angle between the displacement and the direction of the spring, and Ãij = 0 otherwise.The restoring force on each particle is the sum of the forces from all the springs it is connected to, so that i.e. Defining Even though we have presented the argument as if the displacement of each particle is one scalar variable, it is easy to extend this to particles in a d-dimensional system, with N displacement variables for N/d particles.We have implicitly assumed that the particles are frictionless spheres, so that torque balance is trivially satisfied.
The matrix A is a rectangular matrix, since the number of springs is greater than N.As one approaches Point J, the number of contact forces decreases, being equal to N at the transition.
In the random matrix approach to this problem, we assume that all the entries in the matrix A are independent Gaussian random variables, drawn from a distribution with zero mean and (with a suitable rescaling) unit variance.This is the Laguerre random matrix ensemble.It can be shown [24] that the reduced probability density for the eigenvalues ω where the matrix A is a M × N dimensional matrix, i.e. there are M inter-particle contacts in the system, and without loss of generality we have chosen all the ω i 's to be greater than zero.
Rewriting the probability density in terms of wherever ρ(λ) = 0.Here ρ(λ) is the density of eigenvalues, normalized to ∞ 0 ρ(λ)dλ = 1, and the P denotes the principal value of the integral.Symmetrizing ρ(λ) by defining ρ(λ < 0) = ρ(−λ), the function of the complex variable λ is analytic everywhere except that it has branch cuts on the real line over intervals where ρ(λ) = 0, where it is equal to Furthermore, F (λ → 0) is finite, and F (λ → ∞) → 2/λ, because the symmetrized extension of ρ(λ) integrates to 2. This has the solution The density of states is then where we have removed the extension of ρ(λ) to λ < 0, so that has the same form, but with a and b rescaled, which is not significant since Eq.( 5) was already obtained after rescaling.
When M/N > 1, there is a broad peak in D(ω), with a gap in the spectrum near ω = 0.The peak is not symmetric, falling off much more sharply on the small ω side than on the large ω side.In the middle, the peak slopes downwards as ω is increased.As M/N is reduced, the gap shrinks while the width of the peak remains constant.
One can compare these analytical predictions with numerical results.Simulations of two-dimensional frictionless soft spheres are performed, allowing the spheres to equilibrate from an initial random configuration [25].As with the analytical prediction, there is a broad peak, that falls off more sharply at small ω than at large ω.The density of states D(ω = 0) = 0 except at the transition.However, the numerical data does not show the gap in the spectrum near ω = 0 predicted by random matrix theory.The numerical data also has a pronounced boson peak at a non-zero value of ω, and a cusp in D(ω) at the origin at the transition.None of these is consistent with the prediction from random matrix theory.We return to this point in Section IV.

III. CORRELATIONS
In this section we subject the predictions of the random matrix model to more stringent and appropriate tests.We have seen in the previous section the mean density of states of the random matrix model does not exactly match the density of states of the dynamical numerical simulation of a granular pack.However that is not the appropriate test of a random matrix model.In all of its successful applications, what random matrix theory is able to predict correctly is not the mean density of states, but rather statistical features like the density of states correlation and the level spacing distribution, that are computed after the spectrum has been smoothed by a procedure called unfolding (explained below) that makes the mean density of states uniform.
Theoretically one can understand this as follows.It can be shown [26] that, for an ensemble of random matrices, the mean density of states can be changed at will by varying the assumed distribution of the matrix elements, but that the correlations of the unfolded spectrum retain a universal form that depends only on the symmetries of the ensemble considered (such as whether the matrices are real or complex or quaternion real).From this point of view the mean density of states we obtained in the previous section is simply an artifact of the Gaussian distribution we chose for our random matrix elements, but the unfolded density of states correlations and the level spacing distribution are universal and should match the numerical data for jammed granular matter if the model FIG. 1.The red and blue histograms show the consecutive level spacing distribution for the numerically computed vibrational spectra of jammed granular material.An ensemble of one thousand realizations of the jammed state was used.The red histogram bins the eleven consecutive level spacings between the frequencies ω5 through ω16 for each realization; the blue histogram eleven consecutive spacings between frequencies ω400 through ω411.Each spacing ωi+1 − ωi is normalized by ωi+1 − ωi , where the average is taken over the one thousand realizations.The black curve corresponds to the Wigner surmise for the level spacing distribution of the Gaussian orthogonal ensemble (which is indistinguishable from the Laguerre ensemble at this level of resolution).The close agreement between the two histograms and the solid black curve are consistent with the predictions of our random matrix model of the jammed state of granular matter.

is applicable.
The key feature of the random matrix spectrum is that it is rigid (i.e.highly correlated).The rigidity of the spectrum is revealed at small energy scales by the distribution of consecutive level spacings.The longer range rigidity can be demonstrated by the autocorrelation of the density of states or by specific statistical measures such as the number statistic and the spectral rigidity [18].
Insight into the strong correlations between the eigenvalues implied by the Laguerre ensemble distribution Eq.( 5) is provided by the following plasma analogy.We focus on the case M = N + 1 since we are interested in the spectrum for point J.If we rewrite rewrite the factor |ω 2 j − ω 2 i | in Eq.( 5) as exp(ln |ω j − ω i | + ln |ω j + ω i |) and the factor ω i as exp(ln ω i ), we can interpret Eq. ( 5) as the partition function of a classical plasma of N particles located on the positive ω axis at the locations ω 1 , . . ., ω N with logarithmic interactions between the particles as well as logarithmic interactions between each particle and image particles at locations −ω 1 , . . ., −ω N .In the plasma analogy the particles are also confined near the origin by a quadratic potential and are constrained to remain on the positive ω axis by a hard wall at the origin.
We turn now to the distribution of spacings between consecutive levels.Intuitively one might expect that the spacing between consecutive levels at low frequency would be different from that between two levels at high frequency.This is because, in terms of the plasma analogy, in the first instance the two interacting levels are near the hard wall, whereas in the latter instance they are deep in the interior of the plasma.However, we show in Appendix A that the consecutive level spacing distribution is not noticeably different for high and low frequencies, and neither of these is noticeably different from the distribution for the GOE.
In Fig. 1 we compare the predictions of the Laguerre model to the numerically computed jammed granular spectrum.The red histogram bins the eleven consecutive level spacings between the frequencies ω 5 through ω 16 for each realization; each spacing is normalized by ω i+1 −ω i , where the average is taken over the one thousand realizations of the jammed granular pack.The blue histogram is the same but for the eleven level spacings between the frequencies ω 400 through ω 411 .The two distributions are seen to be indistinguishable and to be in good agreement with the approximate analytic formula for the Laguerre ensemble (solid black curve) which is itself indistinguishable from the prediction of the GOE at the resolution of the figure.
Fig. 1 shows that the numerical level spacing distribution is consistent with the prediction of the Laguerre random matrix model and hence is a validation of that model.However, the level spacing distribution is not able to distinguish between the Laguerre model and the GOE, and is therefore a less sharp test of the Laguerre model than the density of states autocorrelation to which we now turn.
In order to calculate these correlations it is convenient to rewrite the distribution in Eq. ( 5) for the case M = N + 1 in the form where x i = ω 2 i .The one-point and two-point correlation functions are defined as The plasma analogy shows that calculation of the correlation functions in Eq. ( 14) and Eq. ( 15) is a formidable problem in classical statistical mechanics.Nonetheless it has been exactly done by Nagao and Slevin [19] by rewriting Eq. ( 13) in the form of a quaternion determinant and performing the integrals by a generalization of a theorem of Dyson [27] on integration over quaternion determinants.Before we give those results we first describe the unfolding procedure.R 1 (x) is evidently the density of states, and we now define where ξ(x) is the cumulative density of states.The unfolded two point correlation function is then defined as It is easy to see that if L 1 (ξ) is similarly defined then L 1 (ξ) = 1 showing that after unfolding the density of states is uniform with a mean level spacing of unity.The exact expression for L 2 is rather lengthy and is relegated to Appendix B.
In Figure 2 we plot 1 − L 2 (ξ, 0) as a function of ξ for the Laugerre ensemble.The corresponding plot for the Gaussian Orthogonal ensemble is also shown.When ξ is large, the two curves approach each other.Indeed, the analytical expression for 1 − L 2 (ξ) for ξ 1 in the Laguerre ensemble can be verified to be which coincides with the form for the same quantity in the GOE.Intuitively, the reason for this coincidence can be understood in terms of the plasma analogy.If the particles are deep in the interior of the plasma the edge effects produced by the image charges are screened and the Laguerre plasma becomes indistinguishable from the GOE plasma.
Although the two curves coincide in the asymptotic limit ξ 1, Fig 2 shows that there is a range of ξ values where the predictions of the Laguerre ensemble differ significantly from the GOE.Hence comparison to the correlation function for the numerical data for the jammed granular spectrum provides a stringent test that is able to distinguish between the Laguerre ensemble and the GOE.
The numerical data are analyzed as follows.The vibrational frequencies obtained in the numerical simulations from all the 1000 realizations of the jammed state are merged together, and bins are constructed with 200 eigenvalues in each, i.e. there is an average of 0.2 eigenvalues per realization of the jammed state in each bin.Next, we calculate where n i is the number of vibrational frequencies in the i th bin in any given realization, and the average is over the realizations of the jammed state.The histogram of the values obtained for this discretized correlation function are compared with the analytical prediction from the Laugerre ensemble and the GOE, and as seen in Fig. 2, the Laguerre ensemble fits the data very well within the error bars (while the GOE does not).

IV. RANDOM LATTICE MODEL
As discussed earlier, the extent to which the density D(ω) of vibrational frequencies for jammed granular materials agrees with the predictions of random matrix theory is not a good test of the applicability of random matrix theory to these materials, because the distribution of eigenvalues is a non-universal prediction of random matrix theory: if the random matrices are not assumed to be Gaussian, the density of eigenvalues changes.
Nevertheless, there are qualitative discrepancies between the numerically measured D(ω) and the density of eigenvalues {ω i } obtained from random matrix theory, that are worth trying to address.As seen in Eq.( 11), there is a gap in the spectrum of eigenvalues near ω = 0 above the jamming transition, where M > N in Eq. (11).By contrast, in the numerical simulations, D(ω) is only zero at ω = 0 (except at the jamming transition), and increases linearly for small ω.Second, at the jamming transition, D(ω) has a cusp-like peak at ω = 0, while ran-dom matrix theory predicts a flat D(ω) near ω = 0 at the transition.Finally, the numerical D(ω) has a pronounced boson peak at ω = ω 0 = 0, which is not reproduced by random matrix theory.
Various attempts have been made to construct random matrix models that reproduce the properties of granular systems more closely.Ref. [28] studies several different random matrix ensembles and their effect on the structure of eigenmodes with frequencies in the boson peak.Ref. [29] uses weighted Laplacian dynamical matrices to reproduce an intermediate regime in D(ω) (between the boson peak and the low frequency behavior) and ∼ ω 4 scaling of the density of states in this regime.Ref. [30] uses a combination of a random and a regular matrix for the dynamical matrix, to eliminate the gap near ω = 0. Also, Ref. [31] has studied an abstract model that they argue is in the appropriate universality class.
In Eq.( 4), we have assumed that the entries in the matrix A are all independent random variables drawn from a Gaussian distribution.In reality, since the matrix A is supposed to be a mapping from coordinates to contact forces, and only two particles are associated with a contact, only the entries associated with two particles (with d entries per particle for d-dimensional particles) should be non-zero in any column of A. Thus A should be a sparse matrix.
One could choose the two particles associated with each force randomly, but this would result in the system breaking up into separate clusters, not connected to each other, leading to an overabundance of zero modes.Moreover, the concept of adjacency would not be respected: two randomly chosen particles would be likely to be far apart, and should not have been allowed to share a contact.
Instead of choosing the particles associated with a force randomly, we approximate the system as being equivalent to a triangular lattice (with periodic boundary conditions), but with each particle displaced from the position where it would be in a perfect triangular lattice.This randomizes the orientation of the contacts between particles.
To be specific, particles are arranged in successive horizontal layers, with each particle having contacts with the two particles immediately below it: slightly to the left and slightly to the right.Shifting the numbering in each row by half a lattice spacing relative to its predecessor, the particle (i, j) connects to the particles numbered (i, j − 1) and (i + 1, j − 1) with periodic boundary conditions in both directions.(A particle in the bottom layer, (i, 1), connects with (i−L/2, L) and (i+1−L/2, L) in the topmost layer, where L is the number of layers.)In addition, each particle has a probability of connecting to its neighbor on the right in the same row: (i, j) → (i + 1, j).All contacts are bidirectional, i.e. each particle is connected to two particles in the row above it, two in the row below it, and either zero, one or two adjacent particles in the same row.The spring constant associated with each contact is chosen randomly, and the bond angles are also chosen randomly with the constraint that the bond connecting a particle to its left (right) neighbor in the row below points down and to the left (right), while the bond connecting a particle to its neighbor in the same row on the left (right) is more horizontal than the bond to its neightbor below and to the left (right).This model is similar to the model introduced for free-standing granular piles [32], a vector generalization of the scalar-force "q-Model" used to model such systems [1].This Random Lattice Model (RLM) with 24 × 24 sites was simulated in this manner, and the vibrational frequencies from 100 different realizations of randomness were merged and plotted as a histogram.In twodimensions, the number of coordinate degrees of freedom is N = 2 × 24 × 24 = 1152, which is comparable to the 800 frequencies that were present in each system in the dynamical simulations [25].The ratio of the number of contact forces to the number of coordinate degrees of freedom, which corresponds to M/N, increases from 1 to 1.5 as the probability of establishing contacts within the same layer increases from 0 to 1.
The results are shown in Figure 3.The transition from D(ω = 0) = 0 when M/N = 1 to D(ω = 0) = 0 when M/N = 1, which is seen in random matrix theory and in the dynamical simulations, is also seen in the RLM density of states.But in addition, there is a cusp in D(ω = 0) at the transition point, as is seen in the dynamical simulations [25].The boson peak at ω = 0 seen in the simulations is also reproduced for the lattice model, being most pronounced at the transition point.
Although there is still a gap in the Random Lattice Model frequency spectrum near ω = 0, it is approximately half what one would predict from random matrix theory for the corresponding M/N.Moreover, it is clear that it is a finite size effect: the global translational invariance of the random lattice with periodic boundary conditions ensures that there are two zero modes (seen clearly in the first plot in Figure 3), and the locality of the connections that are made ensures that long wavelength oscillations have low frequencies.This is verified by increasing N and checking that the gap decreases even though M/N is held constant.
We see that the Random Lattice Model has the same change in D(ω = 0) at the jamming transition that is seen in the dynamical numerical simulations and in random matrix theory.However, it also reproduces other qualitative features of the numerical D(ω) that random matrix theory did not.As seen in Figure 4 and in Figure 5, the distribution of spacings between consecutive frequencies is found to be the same as for the random matrix ensemble, consistent with the Wigner surmise, and the correlation function 1 − L 2 (ξ, 0) at M/N = 1 matches that obtained for the Laguerre ensemble.Thus the random lattice model retains the positive features of random matrix theory, while curing its problems.

V. CONCLUSIONS
In this paper, we show that a random matrix approach can be used successfully to calculate the correlations between vibrational frequencies in a granular system near the jamming transition, if the matrix ensemble is chosen correctly.By modifying the random matrices according to physical considerations, a Random Lattice Model is constructed, which retains the correlation functions of random matrix theory and also successfully reproduces all the qualitative features in the density of vibrational frequencies.Such random lattice models may be more broadly applicable to granular materials.itively Wigner's surmise works because at small spacing the distribution is dominated by the interaction between the two consecutive levels; the other levels that are neglected in the analysis only matter at large spacings.In the same spirit we consider the distribution Eq. ( 5) for the case M = N + 1 = 3.We define ∆ = ω 1 − ω 2 as the spacing and ω = 1 2 (ω 1 + ω 2 ) as the mean energy of the two levels, and we take ω 1 > ω 2 .For a fixed ω the level spacing distribution is then given by for 0 < ∆ < 2ω and p = 0 for ∆ > 2ω.Imposing the condition that the distribution p is normalized and rescaling the level spacing so that ∆ = 1 we obtain the analog of the Wigner surmise for the ensemble in Eq. ( 5).Because all the integrals can be worked out in closed form, an explicit but extremely lengthy formula can be given for p(∆, ω) with appropriate normalization and scaling.
For the sake of brevity we omit this formula but show in Fig. 6 a plot of p(∆, ω) for fixed ω = 5.The distribution is essentially indistinguishable from the Wigner surmise for the Gaussian orthogonal ensemble.The deviation |p(∆, 5.0) − p goe (∆)| < 2 × 10 −4 over the range of the plot and hence invisible to the eye at the resolution of the figure.ω is a measure of how close the two consecutive levels are to the edge; to be precise ω 0 dΩR 1 (Ω) is the ordinal number of the pair whose spacing distribution is being considered; here R 1 is the mean density of states corresponding to the distribution in Eq. ( 5).As one might expect the distribution approaches that for the Gaussian orthogonal ensemble as ω increases.More surprisingly and perhaps disappointingly we find that for pairs of levels that are quite near to zero frequency also the Gaussian orthogonal ensemble is a good approximation.This means that in testing whether the spectrum of jammed granular matter is described by random matrix theory we cannot use the level spacing to distinguish between our model and the Gaussian orthogonal ensemble.Since the derivation of the level spacing distribution above is non-rigorous we have tested it by directly simulating our random matrix model.To this end we generated an ensemble of ten thousand M × N random matrices A with M = N + 1 = 101.The matrix elements are drawn from a Gaussian distribution with zero mean and unit variance.We then evaluated the eigenvalues, ω 2 i of A T A. Fig. 6 shows a histogram of the spacing between the fifth and sixth levels, ω 6 − ω 5 , in the different realizations.Those spacings have been scaled to have unit mean.As can be seen from the figure the level spacing distribution derived here as well as the distribution for the Gaussian orthogonal ensemble provide an excellent fit to the numerical data.Using the numerical data we were also able to test that the distribution for higher levels is practically indistinguishable, confirming all the features of the analysis above.
It is possible to derive an exact expression for the level spacing distribution using the technology of quaternion determinants developed by Dyson [27] and Mehta [18] and generalized to the Laguerre ensemble by Nagao and Slevin [19].This analysis would allow a more careful study of the tails of the distribution where it might deviate from the simple approximation derived here.Such an analysis is not needed for the application considered in this paper but is of intrinsic mathematical interest and we will return to it elsewhere.
Appendix B: Correlations of the Laguerre Ensemble Nagao and Slevin [19] have shown that the correlations for the Laguerre ensemble Eq. ( 13) can be computed by rewriting the distribution as a quaternion determinant and performing the required integrals using a powerful generalization of a theorem by Dyson [27].Here we summarize their results, taking the opportunity to correct some typos in their paper, and to present the results in a form that does not require the reader to have familiarity with the specialized language of quaternion determinants or Pfaffians.
We start by defining the unfolding function (B1) Note that the order of the Bessel function in the first term of the integrand above is given incorrectly in Ref. [19].
Next we define and x .
(B3) In terms of these functions one can write down where it is understood that x is short for x(ξ) and x for x(ξ ), the inverse of the function given by Eq. (B1).We also write B5) where θ denotes the unit step function and the variable of integration and the arguments of the integrand in Eq. (B5) are given incorrectly in Ref. [19].
With these definitions established we can now write down the principal result of Nagao and Slevin [19]  Eqs.(B7) and (B8) are the principal results needed for our test of the random matrix model.Fig. 7 shows a plot of 1 − L 2 (ξ, 0) (solid black curve) as well as the corresponding correlation for the GOE (dotted black curve) which is given by the expression in Eq. (18).We have also computed the correlation function for the ensemble of ten thousand spectra drawn from the Laguerre ensemble as described in in Appendix A. This is shown as the blue histogram.Since these spectra correspond to the very ensemble for which the solid curve represents the exact solution the reason for any departure from the analytic result must be attributed to the finite width of the bins used in analyzing the numerical data and in the finite size of the sample of spectra.This is confirmed by the red histogram which is based on a smaller ensemble of one thousand spectra and shows markedly bigger departures from the analytic result especially for larger ξ.Comparing the red histogram in Fig. 7 to the histogram in Fig. 2 we see that the two agree about equally well with the analytic curve showing that the Laguerre ensemble prediction is indeed in excellent agreement with the numerical data for the vibrational spectrum of jammed granular packs.In more detail the histograms are calculated as follows.We generate an ensemble of M ×N random matrices A α where α = 1, . . ., µ and µ is the number of realizations in the ensemble.We then compute the eigenvalues of A T α A α and denote them ω 2 iα where i = 1 . . .N .The range over which the eigenvalues lie is divided into ν bins and the calculated spectra are binned.Denoting by n iα the number of levels in bin i in realization α we then compute n i is the unfolded width of the bin.The density of states correlator between different bins is then estimated by

FIG. 2 .
FIG.2.Correlation function 1 − L2(ξ, 0) from random matrix theory as a function of ξ, out to ξ = 3.0, i.e. out to a distance within which there should, on average, be three eigenvalues.This is compared to a histogram of the vibrational frequencies from the numerical simulations.The analytical results for the Gaussian orthogonal ensemble and the Laugerre ensemble are both shown, and the latter is clearly superior.

FIG. 3 .
FIG.3.Histogram of vibrational frequencies obtained from the random lattice model described in this paper.100 realizations of 24 x 24 random triangular lattices were constructed, with the ratio M/N equal to 1.1 (top), 1.05 (middle) and 1.0 (bottom) respectively.

FIG. 4 .
FIG. 4. Correlation function 1 − L2(ξ, 0) for the Random Lattice Model introduced in this paper and for the Laguerre ensemble, at the transition point (i.e. with square random matrices).Good agreement is seen between the two.The correlation function for the Gaussian Orthogonal ensemble is also shown for comparison.

FIG. 5 .
FIG.5.Level spacings for the Random Lattice Model and a fit to the Wigner surmise for the Gaussian Orthogonal Ensemble.Spacings from the fifth to the fifteenth normal mode frequencies are normalized, as discussed in the paper, and combined to create the histogram.The Wigner surmise fits the distribution very well, but as discussed in the paper, the fit applies equally to the Laguerre ensemble.

FIG. 6 .
FIG. 6. Plot of the Wigner surmise for the consecutive level spacing distribution for the Laguerre ensemble (solid curve).At the resolution of this plot the change in the distribution as a function of the distance of the levels from the zero frequency edge as well as the deviation from the Gaussian orthogonal ensemble are invisible.The specific curve plotted in the figure corresponds to p(∆, ω) with ω = 5.0; see the discussion in Appendix A. This curve is seen to be in excellent agreement with a histogram obtained from a numerical simulation of the Laguerre ensemble.

namely 1 −
FIG.7.Plot of the correlation function for the Laguerre ensemble (solid curve) and the GOE (dotted curve).The blue histogram is the correlation function for an ensemble of 10,000 spectra drawn from the Laguerre ensemble generated as described in Appendix A. The red histogram is the same but for a smaller ensemble of 1000 spectra.Comparison of the two reveals that the departure of the numerical correlation from the exact analytic result is due in part to the finite bind width (which is the same in both case) but in part also due to the finite size of the ensemble used to estimate the correlation.