Vibrational spectrum of Granular packings with random matrices

TL;DR

This paper shows that the vibrational spectrum of granular materials can be modeled using random matrix theory, specifically the Laguerre orthogonal ensemble.

Contribution

The paper introduces a random lattice model that captures the vibrational properties of granular matter and validates the use of the Laguerre orthogonal ensemble.

Findings

The autocorrelation of the vibrational density of states matches the Laguerre orthogonal ensemble predictions.

The Gaussian orthogonal ensemble fails to reproduce the statistical properties of jammed granular matter.

A random lattice model successfully replicates the vibrational features of frictionless granular matter.

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.

Figures

Plot of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\omega )$$\end{document}$ vs $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega ,$$\end{document}$ with $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{

The correlation function $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1 - L_2 (\xi , 0)$$\end{document}$ is computed using vibrational spectra from 1000 realizations of the jammed state. It appears as a staircase in the plot due to the finite widths of the bins used to compute the correlation. Also shown are the analytic results for the Gaussian orthogonal ensemble (dotted blue curve) and the Laguerre ensemble

The red (solid) and blue (dashed) histograms show the consecutive level spacing distribution for the numerically computed vibrational spectra of jammed frictionless 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 $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _5$$\end{document}$ th

Arrangement of particles in the random lattice model. Each particle rests on at least the two nearest neighbors below it, shown as solid lines. In addition, a particle occasionally also has a contact with a next nearest neighbor below it on one side or the other. For clarity, these are shown as dashed lines, even though they are completely equivalent to the solid bonds. Each particle therefore has two, three or four contacts from above. Periodic boundary conditions are imposed on all four boundaries. As shown, the numbering convention results in the particle (*x*, *y*) being shifted horizontal

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

Correlation function $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1 - L_2(\xi , 0)$$\end{document}$ for the random lattice model introduced in this paper (red staircase) and for the Laguerre ensemble at the transition point (solid black curve). Good agreement is seen between the two. The correlation function for the Gaussian orthogonal ensemble is also shown for comparison (dotted blue curve)

Histogram of level spacings for the random lattice model (red) and a fit to the Wigner surmise for the Gaussian orthogonal ensemble (solid black curve). 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, this applies equally to the Wigner surmise for the Laguerre ensemble

Full text

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, 2] and the continuum limit [3] are difficult to obtain. This is because interparticle contacts are very stiff: a slight compression of two particles that are in 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; [4] and that there are frictional forces between particles, [5] 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. [6] Experimentally, one observes distinctive phenomena such as force chains and stability against mechanical collapse in very sparse static packings, non-Maxwellian velocity distributions [7–9] and inelastic collapse in dilute granular gases, [10, 11] and shear thinning and shear thickening. [12]

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. [13] 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. [14] When the transition occurs at zero temperature and zero shear stress as the density is varied, the transition point is called ‘Point J,’ [15] and is characterized by diverging length scales [16, 17] suggestive of a second-order phase transition. At the same time, other properties of the system change discontinuously at Point J [15] 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 nonzero compression) can be analyzed using a linear model, resulting in normal modes. Extensive numerical simulations [15] on systems at zero temperature and zero shear stress show that the density of states \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\omega )$$\end{document} as a function of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega $$\end{document} approaches zero linearly as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega \rightarrow 0$$\end{document} if the particle density is greater than the critical density. As the particle density is reduced, the slope of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\omega )$$\end{document} at the origin becomes steeper, until at Point J, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\omega \rightarrow 0) \ne 0.$$\end{document}

In the linearized analysis 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 random matrices, and see how the resultant density of states evolves near the transition. This has been done [18, 19] using the Laguerre ensemble instead of the Gaussian orthogonal ensemble (GOE) for the random matrices, in accordance with—as discussed in the next section of this paper—the symmetries of the problem. The model yields a broad peak in the density of states that reaches \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega = 0$$\end{document} as the transition is approached. However, the model also predicts a gap in the density of states near \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega = 0$$\end{document} above the transition, which does not match the numerical results. This is not surprising, because the density of states predicted by random matrix theory is well known to suffer from non-universal effects [20]. For instance, the density of states can be changed at will by varying the assumed distribution of the matrix elements. [21] Instead, the correlations in the density of states and the distribution of level spacings are more reliable indicators of the validity of the random matrix approach [20].

The agreement of the level spacing distribution with the GOE result has been observed elsewhere [22, 23] without reference to the Laguerre ensemble. The distribution of level spacings has been studied for the Laguerre ensemble [24], without comparing to the numerical data for bead packs, but it is unfortunately indistinguishable from the results for the GOE. Taken together, Refs. [22–24] show that the level spacing distribution for bead packs is consistent with random matrix theory, without being able to confirm the appropriate ensemble.

To summarize, the density of states for random matrix ensembles is known to be a non-universal feature; not surprisingly, therefore, the density of states for the simplest version of the Laguerre ensemble shows qualitative differences from the density of states of granular vibrational modes observed in numerical simulations. On the other hand, the level spacing distribution is found to be too universal: there is agreement between the Laguerre ensemble and numerical simulations, but there is equally good agreement between the GOE and numerical simulations. Thus, neither of these can be considered as conclusive evidence for the Laguerre ensemble.

In this paper, we therefore turn to the correlations in the density of states predicted by random matrix theory. We consider the correlation function for the Laguerre ensemble, which differs from that for the GOE near the low-frequency edge of the allowed range of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega $$\end{document} [25]. By comparing the correlations in the numerically computed vibrational spectrum of frictionless 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.

In addition, we 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\omega )$$\end{document} 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. Models that are essentially mean field versions of our model have been studied earlier [26, 27], and differ from our results mainly in the fact that the low-frequency \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\omega ) \sim \omega ^{d-1}$$\end{document} behavior for a d-dimensional system is not seen. The random lattice model is a variation of a model used earlier to study the static properties of free-standing granular piles (with vector forces) [28]; the fact that such similar models yield results in agreement with different experiments increases the credibility of the model.

Other authors have studied variations of random matrix theory: Ref [29] is an early random matrix model for the related problem of the vibrational spectrum of glasses close to the glass transition. Reference [30] studies several different random matrix ensembles and their effect on the structure of eigenmodes with frequencies in the boson peak. Reference [31] uses weighted Laplacian dynamical matrices to reproduce an intermediate regime in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\omega )$$\end{document} (between the boson peak and the low-frequency behavior) and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim \omega ^4$$\end{document} scaling of the density of states in this regime. Reference [32] uses a combination of a random and a regular matrix for the dynamical matrix, to eliminate the gap near \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega = 0.$$\end{document} Also, Ref. [33] has studied an abstract model that they argue is in the appropriate universality class.

The rest of this paper is organized as follows. Section 2 reviews why the symmetry properties of the systems in which we are interested are in the Laguerre universality class rather than the GOE. Section 3 has the first main result of this paper: that the analytically calculated autocorrelation function for the density of energy levels of the Laguerre orthogonal ensemble agrees with the correlation function from dynamical simulations on bead packs, in the vicinity of Point J, but that this is not true for the correlation function for the Gaussian orthogonal ensemble. Section 4 has the second main result of this paper: that the random lattice model, which retains some spatial correlations discarded by the simple random matrix model, is numerically found to have a similar behavior for the density of energy levels as seen in the molecular dynamics simulations, both at and away from Point J. The energy level correlation function and the level spacings for the random lattice model are also seen to agree with the molecular dynamics results (and the Laguerre ensemble results), in accordance with the expectation that such correlations should only depend on the symmetries of the problem instead of the model details. Thus, the random lattice model retains the points of agreement between the Laguerre ensemble and molecular dynamics, and also successfully reproduces the density of states.

Laguerre ensembles

We follow the approach of Ref. [18, 19] here. Within linear response, if the particles in a frictionless granular assembly are displaced slightly from their resting positions, their accelerations are of the form \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ddot{x} = -K x,$$\end{document} where x is a N component column vector with (a rescaled version of) the displacements of all the particles and K is a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\times N$$\end{document} matrix called the dynamical matrix [34]. If d is the dimensionality of the granular assembly, then each particle has d components to its displacement, and N is therefore equal to d times the number of particles.

The crucial observation [35, 36] 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. Thus, we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \Delta \ell _\alpha = {\tilde{A}}_{\alpha j} x_j \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta \ell _\alpha $$\end{document} is the amount by which the spring labeled \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} is compressed. Here, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_j$$\end{document} is one component of the displacement of one of the particles in the assembly. If the particle is one of the two that is involved in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} ’th contact, then \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{A}}_{\alpha j}$$\end{document} is simply \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\cos \theta _{\alpha j},$$\end{document} where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _{\alpha j}$$\end{document} is the angle between the direction of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} ’th spring and the direction of the j’th displacement, i.e., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\cos \theta _{\alpha j} x_j$$\end{document} is the compression of the contact caused by the displacement. If \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_j$$\end{document} is not associated with either of the particles involved in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} ’th contact, then \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{A}_{\alpha j} = 0.$$\end{document} Note that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{A}}$$\end{document} is a rectangular matrix: if there are M interparticle contacts and N particle displacements, then \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{A}}$$\end{document} is an \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M\times N$$\end{document} matrix. As one approaches Point J, the number of contact forces decreases, being equal to N at the transition.

The spring exerts a restoring force that is proportional to this compression; the spring constant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_\alpha $$\end{document} can be different for each spring. The restoring force on each particle is the sum of the forces from all the springs it is connected to. Therefore,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} m_j \ddot{x}_j = - {\tilde{A}}^T_{j\alpha } k_\alpha \Delta \ell _\alpha , \end{aligned}$$\end{document}

i.e., the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} ’th contact force acts as a restoring force on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_j$$\end{document} only if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_j$$\end{document} is associated with one of the two particles involved in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} ’th contact, in which case \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A^T_{j\alpha }$$\end{document} is the same projection factor \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\cos \theta _{\alpha j}$$\end{document} that we had in the dependence of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_\alpha $$\end{document} on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_j.$$\end{document} Putting this equation together with Eq. (1), we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} m_j \ddot{x}_j = - {\tilde{A}}^T_{j\alpha } k_\alpha {\tilde{A}}_{\alpha i} x_i. \end{aligned}$$\end{document}

Defining \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_{\alpha j} = \sqrt{k}_\alpha {\tilde{A}}_{\alpha j}/\sqrt{m}_j,$$\end{document} this is equivalent to

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \sqrt{m_j}\ddot{x}_j = - A^T_{j\alpha } A_{\alpha i} \sqrt{m_i} x_i \end{aligned}$$\end{document}

We have implicitly assumed that the particles are frictionless spheres, so that torque balance is trivially satisfied. Absorbing a factor of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{m_j}$$\end{document} in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_j$$\end{document} for each j,  we finally have the equation

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \ddot{x}_j = - A^T_{j\alpha } A_{\alpha i} x_i , \end{aligned}$$\end{document}

i.e., the dynamical matrix K is a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\times N$$\end{document} matrix that is equal to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A^T A.$$\end{document}

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 (also known as the Wishart ensemble in the mathematical literature), and is different from assuming that the elements of K are independent Gaussian random variables (with the constraint that K is a symmetric matrix). The approach to Point J is modeled by adjusting the ratio M/N to approach 1 from above.

The random matrix approach ignores the fact that each contact force acts between only two particles, i.e., that each row in the matrix A has only 2d entries. It also ignores the fact that the network of contacts is created by particles rearranging themselves, which would be expected to result in correlations between the contact forces. To address these limitations, one would need to construct a model which includes spatial information about the contacts (or perform a full molecular dynamics simulation), discussed in Sect. 4 of this paper. When this is done, the non-universal results from random matrix theory are affected without changing the universal results.

It can be shown [37] that, for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M \ge N,$$\end{document} the density of states for the eigenfrequencies \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _1, \omega _2\ldots \omega _N$$\end{document} is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \rho (\omega ){} & {} = \frac{1}{N \pi } \frac{\sqrt{(N b^2 - \omega ^2)(\omega ^2 - N a^2)}}{\omega } \nonumber \\{} & {} \quad a \sqrt{N}< \omega < b\sqrt{N} \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = \sqrt{M/N}-1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b = \sqrt{M/N} + 1.$$\end{document} When \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M/N > 1,$$\end{document} there is a broad peak in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\omega ),$$\end{document} with a gap in the spectrum near \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega = 0.$$\end{document} The peak is not symmetric, falling off much more sharply on the small \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega $$\end{document} side than on the large \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega $$\end{document} side. In the middle, the peak slopes downwards as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega $$\end{document} is increased. As M/N is reduced, the gap shrinks while the width of the peak remains constant. When \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M/N = 1,$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\omega ) = \sqrt{4 N - \omega ^2}/(N\pi )$$\end{document} which matches the Wigner semicircle law for the Gaussian orthogonal ensemble, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\omega =0)\ne 0.$$\end{document} Fig. 1. Plot of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\omega )$$\end{document} vs \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega ,$$\end{document} with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega $$\end{document} scaled by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{N},$$\end{document} for (red, solid) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M/N = 1.0$$\end{document} and (blue, dashed) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M/N = 1.2.$$\end{document} For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M/N = 1.2,$$\end{document} the gap at small \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega $$\end{document} and the asymmetry of the peak are clearly visible. As M/N is increased, the asymmetry in the shape persists, the gap at small \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega $$\end{document} broadens, and the range of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega /\sqrt{N}$$\end{document} for which \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\omega ) \ne 0$$\end{document} remains equal to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b-a = 2$$\end{document}

One can compare these analytical predictions with numerical results, with data from the O’Hern group [38]. The details of how the numerical simulations were conducted are given in Ref. [38], but to summarize: approximately 1000 frictionless bidisperse disks (half large and half small) with a diameter ratio of 1.4 were placed randomly in a square cell with periodic boundary conditions, and allowed to relax to equilibrium, at which point the frequencies of the normal modes were measured. The interaction between the disks was harmonic (only under compression), but the spring constant was different for large–large, large–small and small–small contacts. As with the analytical prediction, there is a broad peak that falls off more sharply at small \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega $$\end{document} than at large \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega .$$\end{document} The density of states \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\omega =0) = 0$$\end{document} except at the transition. However, the numerical data does not show the gap in the spectrum near \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega = 0$$\end{document} predicted by random matrix theory. The numerical data also has a pronounced boson peak at a nonzero value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega ,$$\end{document} and a cusp in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\omega )$$\end{document} at the origin at the transition [15]. These features do not agree with the prediction from random matrix theory, but as already discussed, random matrix theory is not expected to yield the (non-universal) density of states correctly. We return to this point in Sect. 4, on the random lattice model.

Correlations

In this section, we subject the predictions of the random matrix model to more stringent and appropriate tests, by studying the (universal) correlations. 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 [20].

Insight into the strong correlations between the eigenvalues implied by the Laguerre ensemble distribution is provided by the following plasma analogy. We focus on the case \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M = N+1$$\end{document} since we are interested in the spectrum for point J. For the Laguerre ensemble, the rigidity of the spectrum arises from the factors of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|\omega _j^2 - \omega _i^2|$$\end{document} in the joint probability distribution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p (\omega _1, \omega _2, \ldots , \omega _N)$$\end{document} for the eigenvalues given by Eq. (A2) If we rewrite this as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\exp (\ln |\omega _j - \omega _i | + \ln |\omega _j + \omega _i |),$$\end{document} we can interpret the probability distribution as the partition function of a classical plasma of N particles located on the positive \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega $$\end{document} axis at the locations \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _1, \ldots , \omega _N$$\end{document} with logarithmic interactions between the particles as well and between each particle and image particles at locations \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-\omega _1, \ldots , -\omega _N$$\end{document} . The particles are also confined near the origin by a quadratic potential and are constrained to remain on the positive \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega $$\end{document} axis by a hard wall at the origin. The plasma analogy is another argument to only consider correlations for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M\approx N.$$\end{document} If \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M/N > 1$$\end{document} (in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\rightarrow \infty $$\end{document} limit), the gap in the Laguerre spectrum \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\omega )$$\end{document} near \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega = 0$$\end{document} means that the particles and their images are well separated. The correlations between energy levels in the Laguerre ensemble will then be indistinguishable from those for the GOE. Even for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M\approx N,$$\end{document} the difference between the two ensembles should be greatest near \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega = 0.$$\end{document}

In terms of the variables \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_i = \omega _i^2$$\end{document} , the one-point and two-point correlation functions are defined as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} R_1 (x) = N \int _0^\infty d x_2 \ldots \int _0^\infty d x_N P(x, x_2, \ldots , x_N)\nonumber \\ \end{aligned}$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}{} & {} R_2 (x, y) = N(N-1) \int _0^\infty d x_3 \ldots \nonumber \\{} & {} \quad \int _0^\infty d x_N P(x, y, x_3,\ldots , x_N) \end{aligned}$$\end{document}

Here, P denotes the joint probability distribution of the squared frequencies \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_i = \omega _i^2$$\end{document} which can easily be deduced from the joint probability distribution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p ( \omega _1, \omega _2, \ldots , \omega _N)$$\end{document} given in Eq. (A2). The plasma analogy shows that calculation of the correlation functions in Eqs. (7) and (8) is a formidable problem in classical statistical mechanics. Nonetheless, it has been exactly done by Nagao and Slevin [25] by rewriting the probability density in the form of a quaternion determinant and performing the integrals by a generalization of a theorem of Dyson [39] on integration over quaternion determinants. Before we give those results, we first describe the unfolding procedure.

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_1(x)$$\end{document} is evidently the density of states, and we now define

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \xi (x) = \int _0^x d x^\prime \; R_1 (x^\prime ) \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi (x)$$\end{document} is the cumulative density of states. The unfolded two point correlation function is then defined as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} L_2 (\xi _1, \xi _2) = \frac{1}{R_1[ x (\xi _1) ]} \frac{1}{R_1 [ x(\xi _2) ]} R_2 [ x (\xi _1), x(\xi _2)]\nonumber \\ \end{aligned}$$\end{document}

The exact expression for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_2$$\end{document} is rather lengthy and is given in the Appendix.

In Fig. 2, we plot \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1 - L_2(\xi , 0)$$\end{document} as a function of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi $$\end{document} for the Laguerre ensemble. The corresponding plot for the Gaussian orthogonal ensemble is also shown. When \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi $$\end{document} is large, the two curves approach each other. Indeed, the analytical expression for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1 - L_2(\xi )$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi \gg 1$$\end{document} in the Laguerre ensemble can be verified to be

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}{} & {} 1 - L_2 (\xi , 0)\nonumber \\ {}{} & {} = \frac{\sin ^2 \pi \xi }{ ( \pi \xi )^2} + \left[ \frac{1}{2} - \int _0^\xi d y \frac{ \sin \pi y}{\pi y} \right] \left[ \frac{ \cos \pi \xi }{\xi } - \frac{1}{\pi \xi ^2} \sin \pi \xi \right] \nonumber \\ \end{aligned}$$\end{document}

which coincides with the form for the same quantity in the GOE, as expected since the effect of the image charges in the plasma analogy for the Laguerre ensemble should be small when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi $$\end{document} is large. However, although the two curves coincide in the asymptotic limit \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi \gg 1$$\end{document} , Fig. 2 shows that there is a range of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi $$\end{document} 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 spectrum of the jammed bead pack provides a stringent test that is able to distinguish between the Laguerre ensemble and the GOE.Fig. 2. The correlation function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1 - L_2 (\xi , 0)$$\end{document} is computed using vibrational spectra from 1000 realizations of the jammed state. It appears as a staircase in the plot due to the finite widths of the bins used to compute the correlation. Also shown are the analytic results for the Gaussian orthogonal ensemble (dotted blue curve) and the Laguerre ensemble (solid black curve). The agreement of the vibrational data to the Laguerre ensemble is clearly superior. The vibrational spectrum data are based on simulations by Kyle Vander Werf and Corey O’Hern described in Sect. 2

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

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} 1 - \frac{1}{(0.2)^2} [\langle n_0 n_i\rangle - \langle n_0 \rangle \delta _{i0}] \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_i$$\end{document} is the number of vibrational frequencies in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i^{\textrm{th}}$$\end{document} 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 is compared with the analytical prediction from the Laguerre 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).

For completeness, we also compute the distribution of level spacings from the numerically computed vibrational spectra for bead packs [38]. This is shown in Fig. 3. Although the plasma analogy would suggest that the level spacing distribution could be different in the vicinity of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega = 0$$\end{document} and for large \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega ,$$\end{document} this is not found to be the case; both distributions are found to be indistinguishable from the Wigner surmise, which also matches—within our level of resolution— the distributions for the Laguerre ensemble [24] and the GOE. Thus, although the level spacing distribution is consistent with the Laguerre ensemble, it does not discriminate between the Laguerre and Gaussian orthogonal ensembles.Fig. 3. The red (solid) and blue (dashed) histograms show the consecutive level spacing distribution for the numerically computed vibrational spectra of jammed frictionless 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _5$$\end{document} through \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _{16}$$\end{document} for each realization; the blue histogram eleven consecutive spacings between frequencies \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _{400}$$\end{document} through \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _{411}$$\end{document} . Each spacing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _{i+1} - \omega _i$$\end{document} is normalized by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle \omega _{i+1} - \omega _i\rangle ,$$\end{document} 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 Wigner surmise for the Laguerre ensemble at this level of resolution, as discussed in Appendix Appendix B). 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 frictionless granular matter

Random Lattice model

As discussed earlier, the extent to which the density \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\omega )$$\end{document} of vibrational frequencies for jammed frictionless 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. Nevertheless, there are qualitative discrepancies between the numerically measured \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\omega )$$\end{document} [38] and the density of eigenvalues \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{\omega _i\}$$\end{document} obtained from random matrix theory that are worth trying to address: the low-frequency gap, the absence of a cusp at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega = 0$$\end{document} at the transition and the absence of a boson peak.

In Eq. (5), 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 nonzero in any column of A. Thus, A should be a sparse matrix.

In mean field theory, one would choose the two particles associated with each force randomly. This would result in the system breaking up into separate clusters, not connected to each other, leading to an overabundance of zero modes unless one only retains the giant connected cluster. 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 (Fig. 4).Fig. 4. Arrangement of particles in the random lattice model. Each particle rests on at least the two nearest neighbors below it, shown as solid lines. In addition, a particle occasionally also has a contact with a next nearest neighbor below it on one side or the other. For clarity, these are shown as dashed lines, even though they are completely equivalent to the solid bonds. Each particle therefore has two, three or four contacts from above. Periodic boundary conditions are imposed on all four boundaries. As shown, the numbering convention results in the particle (x, y) being shifted horizontally by half a lattice spacing compared to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(x, y - 1).$$\end{document} The bond angles are chosen randomly, and the vibrational spectrum calculated

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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(i, j - 1)$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(i + 1, j - 1)$$\end{document} with periodic boundary conditions in both directions. (A particle in the bottom layer, (i, 1),  connects with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(i - L/2, L)$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(i + 1 - L/2, L)$$\end{document} in the topmost layer, where L is the number of layers.) In addition, each particle has a probability p of connecting to an additional particle in the row below it, i.e., a probability of 0.5p of connecting to the particle at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(i - 1, j - 1)$$\end{document} and a probability of 0.5p of connecting to the particle at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(i + 2, j - 1),$$\end{document} again with periodic boundary conditions. All contacts are bidirectional, so that each particle is connected to two to four particles in the row above it. The spring constant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa $$\end{document} associated with each contact is an independent random variable, with probability density

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \rho (\kappa )d\kappa = 3 \kappa ^2 \exp [-\kappa ^3] d\kappa . \end{aligned}$$\end{document}

This is what one would expect with exponentially distributed contact forces. Since the contact force f at a Hertzian contact is related to the compression \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta x$$\end{document} of the contact as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f\sim (\delta x)^{3/2},$$\end{document} and the linear stiffness if the compression \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta x$$\end{document} is changed by a small amount is found by differentiating this expression to be \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa \sim (\delta x)^{1/2},$$\end{document} therefore \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f \sim \kappa ^3.$$\end{document} We choose units in which \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f = \kappa ^3.$$\end{document} Thus, an exponential probability density for f implies

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \rho (f) df = \exp [-f] df = \exp [-\kappa ^3] (3 \kappa ^2 d\kappa ) \end{aligned}$$\end{document}

The bond connecting a particle to its left (right) neighbor in the row below points down and to the left (right), and the bond angles are uniform random variables in the third (fourth) quadrant. If the particle at (i, j) is also connected to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(i-1, j - 1),$$\end{document} this bond angle is also chosen to be a uniform random variable in the third quadrant, with the constraint that it must be more horizontal than the bond between (i, j) and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(i, j - 1).$$\end{document} (There is a similar condition if the additional contact is to the right, to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(i + 2, j - 1).$$\end{document} ) This model is similar to a model introduced for free-standing granular piles [28], which is a vector generalization of the scalar-force ‘q-Model’ used to model such systems [1, 2]. We expect that our results will not depend on the details of the model, as long as the spring constants and bond angles are random, the locality of interactions is respected, and the number of contacts made by an individual particle can vary in order to approach the jamming transition.

Having constructed the lattice, it is easy to obtain the matrix A, and therefore the eigenvalues of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A^T A.$$\end{document} One difference from the Laguerre ensemble is that each contact is associated with two particles, and therefore four displacement variables, i.e., each row of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_{\alpha i}$$\end{document} has exactly four nonzero entries. Another difference is that the translational invariance of the system is respected by the matrix A; the consequence of this is discussed later in this section.

Unlike the dynamical simulations, we compare with [38], the random lattice model uses a Hertzian contact force law rather than a harmonic law (with a distribution of spring constants). It also has the property that every lattice site is connected to at least two above it and two below, whereas granular assemblies will have some particles with less than four contacts. Most importantly, the random lattice model does not enforce force balance at each lattice site. The objective with the random lattice model is to show that randomness and locality are sufficient to yield the qualitative aspects of the density of states while preserving the level spacing distribution and the spectral correlations of Random Matrix Theory with the Laguerre ensemble. An improved version of the model, which would be a more complicated version of the model in Ref. [28]—because we do not have gravity to break the symmetry of the system, finding the equilibrium state is more complicated—can be studied if one is interested in how the shape and location of an eigenmode is related to force chains in the underlying stresses. We leave this for future work.

The random lattice model (RLM) with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$128 \times 128$$\end{document} sites was simulated in this manner, and the vibrational frequencies from 100 different realizations of randomness were merged and plotted as a histogram. The ratio of the number of contact forces to the number of coordinate degrees of freedom, which corresponds to M/N,  increases with the probability of a particle establishing extra bonds: from 1 at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p=0$$\end{document} to 1.5 at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p=1.$$\end{document} Fig. 5. Histogram of vibrational frequencies obtained from the random lattice model described in this paper. One hundred realizations of 128x 128 random triangular lattices were constructed, with the ratio M/N equal to 1.1 (top), 1.05 (middle) and 1.0 (bottom) respectively

The results are shown in Fig. 5. When \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M/N = 1,$$\end{document} the density of states \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\omega )$$\end{document} is nonzero at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega = 0,$$\end{document} in agreement with the dynamical simulations [38] and random matrix theory. But in addition, there is a cusp in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\omega = 0)$$\end{document} at the transition point, as seen in the dynamical simulations but not in random matrix theory. (This cusp is also seen in some earlier mean field models [18, 27] which do not have local correlations; the random lattice model also has the advantage of being very similar to a model used successfully to explain force chains in static free-standing sandpiles.) The boson peak at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega \ne 0$$\end{document} seen in the simulations is also reproduced for the lattice model, being most pronounced at the transition point.

Moving away from the transition point at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M/N = 1,$$\end{document} there is a small gap in the spectrum of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\omega )$$\end{document} near \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega = 0.$$\end{document} But this is a finite size effect; we have verified by increasing N while keeping M/N constant that the gap decreases. This is inevitable: the global translational invariance of the random lattice model with periodic boundary conditions ensures that there are two zero modes, and the locality of the connections that are made ensures that long wavelength oscillations have low frequencies. Thus, when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\rightarrow \infty $$\end{document} at constant M/N,  the density of states must scale as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega ^{d-1}$$\end{document} for small \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega .$$\end{document} This is true in the dynamical simulations [38], and would be true for any local model, not just the random lattice model, but it is not true for random matrix theory as seen in Eq. (6) with fixed M/N and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\rightarrow \infty .$$\end{document} Fig. 6. Correlation function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1 - L_2(\xi , 0)$$\end{document} for the random lattice model introduced in this paper (red staircase) and for the Laguerre ensemble at the transition point (solid black curve). Good agreement is seen between the two. The correlation function for the Gaussian orthogonal ensemble is also shown for comparison (dotted blue curve)

We see that the random lattice model reproduces the qualitative features of the numerical density of states \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\omega )$$\end{document} —the cusp in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\omega )$$\end{document} near the origin when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M/N = 1,$$\end{document} the boson peak, the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (\omega )\sim \omega ^{d-1}$$\end{document} scaling when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M/N > 1$$\end{document} —that random matrix theory is unable to do. In addition, as seen in Fig. 6 and in Fig. 7, 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1 - L_2(\xi ,0)$$\end{document} at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M/N = 1$$\end{document} matches that obtained for the Laguerre ensemble. (The correlation function is not compared away from the jamming transition because, as discussed earlier, it is indistinguishable from the GOE correlation function there.) This is what one would hope, since the system is in the universality class of the Laguerre ensemble rather than the GOE (since the dynamical matrix has the form \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A^T A$$\end{document} ), and any reasonable model of the system should be in the same universality class. Thus, it is reassuring that the random lattice model retains the positive features of random matrix theory, while curing its problems.Fig. 7. Histogram of level spacings for the random lattice model (red) and a fit to the Wigner surmise for the Gaussian orthogonal ensemble (solid black curve). 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, this applies equally to the Wigner surmise for the Laguerre ensemble

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. The random lattice model closely resembles a vector generalization of the q-model that has previously been used successfully to understand the distribution of stress in static granular matter [1, 2, 28]. That the same model is able to reproduce both the static and vibrational properties of granular matter suggests it may be more broadly applicable to provide a unified understanding of the physics of granular matter.