Vibrational Spectrum of Granular Packings With Random Matrices

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 Coppersmith and the continuum limit Cates 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; MvHreview and that there are frictional forces between particles, Bi 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. Goldhirsch Experimentally, one observes distinctive phenomena such as force chains and stability against mechanical collapse in very sparse static packings, non-Maxwellian velocity distributions vanNoije ; Rouyer ; Reis and inelastic collapse in dilute granular gases, McNamara ; Hopkins and shear thinning and shear thickening in the intermediate regime. Brown

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. LiuNagelCRCBook 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. LiuNagel1998 When the transition occurs at zero temperature and zero shear stress as the density is varied, the transition point is called “Point J”, OHern and is characterized by diverging length scales Wyart ; Ellenbroek suggestive of a second order phase transition. At the same time, other properties of the system change discontinuously at Point J, OHern 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 OHern on systems at zero temperature and zero shear stress show that the density of states $D(\omega)$ as a function of $\omega$ approaches zero linearly as $\omega\rightarrow 0$ if the particle density is greater than the critical density. As the particle density is reduced, the slope of $D(\omega)$ at the origin becomes steeper, until at Point J, $D(\omega\rightarrow 0)\neq 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 random matrices, and see how the resultant density of states evolves near the transition. This has been done, beltukov and yields a broad peak in the density of states that reaches $\omega=0$ as the transition is approached. However, the model also predicts a gap in the density of states near $\omega=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(\omega)$ 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 mehta . 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 mehta .

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 $\omega$ nagao . 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, Silbert ; Nelson 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(\omega)$ 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.

We follow the approach of Ref. beltukov here. Within linear response, if the particles in the granular assembly are displaced slightly from their resting positions, their accelerations are of the form $\ddot{x}=-Kx,$ where $x$ is a $N$ component column vector and $K$ is a $N\times N$ matrix. The crucial observation lubensky ; calladine 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 $\tilde{A}_{ij}=\cos\theta_{ij}$ if the $i$’th spring is connected to the $j$’th particle, with $\theta_{ij}$ being the angle between the displacement and the direction of the spring, and $\tilde{A}_{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 $A_{ij}=\sqrt{k}_{i}\tilde{A}_{ij}/\sqrt{m}_{j},$ this is equivalent to

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 Forrester that the reduced probability density for the eigenvalues $\omega_{1}^{2},\omega_{2}^{2}\ldots\omega^{2}_{N}$ of $A^{T}A$ is of the form

where the matrix $A$ is a $M\times N$ dimensional matrix, i.e. there are $M$ inter-particle contacts in the system, and without loss of generality we have chosen all the $\omega_{i}$’s to be greater than zero.

Rewriting the probability density in terms of $\lambda_{i}=\omega_{i}/\sqrt{N},$ we have

where $M/N=1+f.$ If $M-N$ is $O(N),$ a saddle-point expansion yields

wherever $\rho(\lambda)\neq 0.$ Here $\rho(\lambda)$ is the density of eigenvalues, normalized to $\int_{0}^{\infty}\rho(\lambda)d\lambda=1,$ and the $P$ denotes the principal value of the integral. Symmetrizing $\rho(\lambda)$ by defining $\rho(\lambda<0)=\rho(-\lambda),$ the function

of the complex variable $\lambda$ is analytic everywhere except that it has branch cuts on the real line over intervals where $\rho(\lambda)\neq 0,$ where it is equal to

Furthermore, $F(\lambda\rightarrow 0)$ is finite, and $F(\lambda\rightarrow\infty)\rightarrow 2/\lambda,$ because the symmetrized extension of $\rho(\lambda)$ integrates to 2. This has the solution

with

The density of states is then

where we have removed the extension of $\rho(\lambda)$ to $\lambda<0,$ so that $\int_{0}^{\infty}\rho(\lambda)d\lambda=1.$ The density of states $D(\omega)$ 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(\omega),$ with a gap in the spectrum near $\omega=0.$ The peak is not symmetric, falling off much more sharply on the small $\omega$ side than on the large $\omega$ side. In the middle, the peak slopes downwards as $\omega$ is increased. As $M/N$ is reduced, the gap shrinks while the width of the peak remains constant. When $M/N=1,$ $D(\omega)=\sqrt{\Omega^{2}-\omega^{2}}/\pi$ which matches the Wigner semicircle law for the Gaussian orthogonal ensemble, and $D(\omega=0)\neq 0.$

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 Corey . As with the analytical prediction, there is a broad peak, that falls off more sharply at small $\omega$ than at large $\omega.$ The density of states $D(\omega=0)=0$ except at the transition. However, the numerical data does not show the gap in the spectrum near $\omega=0$ predicted by random matrix theory. The numerical data also has a pronounced boson peak at a non-zero value of $\omega,$ and a cusp in $D(\omega)$ 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.

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 brezin 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 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 mehta .

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 $|\omega_{j}^{2}-\omega_{i}^{2}|$ in Eq.(5) as $\exp(\ln|\omega_{j}-\omega_{i}|+\ln|\omega_{j}+\omega_{i}|)$ and the factor $\omega_{i}$ as $\exp(\ln\omega_{i})$, we can interpret Eq. (5) as the partition function of a classical plasma of N particles located on the positive $\omega$ axis at the locations $\omega_{1},\ldots,\omega_{N}$ with logarithmic interactions between the particles as well as logarithmic interactions between each particle and image particles at locations $-\omega_{1},\ldots,-\omega_{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 $\omega$ 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 $\omega_{5}$ through $\omega_{16}$ for each realization; each spacing is normalized by $\langle\omega_{i+1}-\omega_{i}\rangle,$ 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 $\omega_{400}$ through $\omega_{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}=\omega_{i}^{2}$. The one-point and two-point correlation functions are defined as

and

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 nagao by rewriting Eq. (13) in the form of a quaternion determinant and performing the integrals by a generalization of a theorem of Dyson dyson 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 $\xi(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}(\xi)$ is similarly defined then $L_{1}(\xi)=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}(\xi,0)$ as a function of $\xi$ for the Laugerre ensemble. The corresponding plot for the Gaussian Orthogonal ensemble is also shown. When $\xi$ is large, the two curves approach each other. Indeed, the analytical expression for $1-L_{2}(\xi)$ for $\xi\gg 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 $\xi\gg 1$, Fig 2 shows that there is a range of $\xi$ 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^{{\rm 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).

As discussed earlier, the extent to which the density $D(\omega)$ 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(\omega)$ and the density of eigenvalues $\{\omega_{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 $\omega=0$ above the jamming transition, where $M>N$ in Eq.(11). By contrast, in the numerical simulations, $D(\omega)$ is only zero at $\omega=0$ (except at the jamming transition), and increases linearly for small $\omega.$ Second, at the jamming transition, $D(\omega)$ has a cusp-like peak at $\omega=0,$ while random matrix theory predicts a flat $D(\omega)$ near $\omega=0$ at the transition. Finally, the numerical $D(\omega)$ has a pronounced boson peak at $\omega=\omega_{0}\neq 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. Manning studies several different random matrix ensembles and their effect on the structure of eigenmodes with frequencies in the boson peak. Ref. Middleton uses weighted Laplacian dynamical matrices to reproduce an intermediate regime in $D(\omega)$ (between the boson peak and the low frequency behavior) and $\sim\omega^{4}$ scaling of the density of states in this regime. Ref. Beltukov1 uses a combination of a random and a regular matrix for the dynamical matrix, to eliminate the gap near $\omega=0.$ Also, Ref. Parisi 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)\rightarrow(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 Narayan , a vector generalization of the scalar-force “q-Model” used to model such systems Coppersmith .

This Random Lattice Model (RLM) with $24\times 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 two-dimensions, the number of coordinate degrees of freedom is $N=2\times 24\times 24=1152,$ which is comparable to the 800 frequencies that were present in each system in the dynamical simulations Corey . 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(\omega=0)=0$ when $M/N\neq 1$ to $D(\omega=0)\neq 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(\omega=0)$ at the transition point, as is seen in the dynamical simulations Corey . The boson peak at $\omega\neq 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 $\omega=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(\omega=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(\omega)$ 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}(\xi,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.

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.