Mean-Field Predictions of Scaling Prefactors Match Low-Dimensional Jammed Packings

Introduction – Granular materials exhibit universal properties regardless of the material properties of the individual grains Liu and Nagel (2010); Charbonneau et al. (2012, 2017). The jamming transition is a critical point near which properties such as pressure, packing fraction, or number of excess contacts, among others, scale as power laws. Scaling theory summarizes and condenses these power law relationships, but no first-principles theory of jammed systems at finite dimensions exists. The replica mean-field theory of glasses and jamming has been shown to be exact in the infinite dimensional limit Parisi and Zamponi (2010); Parisi et al. (2020). To do so it relies on the assumption that there are no correlations between neighbors, fundamentally at odds with low-dimensional systems. As such, mean-field predictions should not be expected to hold in low dimensional-jamming, and some results, most notably the packing fraction at jamming, deviate from the mean-field predictions Charbonneau et al. (2012); Parisi et al. (2018). However, despite the fact that low dimensional systems have highly correlated neighbors the scaling relations are precisely the same as those found in infinite dimensions O’Hern et al. (2003, 2002); Goodrich et al. (2016). Many other results predicted by the mean field have also been observed in low dimensional jamming, suggesting that they may be provable without the mean field approximation Charbonneau et al. (2012, 2016); Berthier et al. (2019); Charbonneau et al. (2017); Dennis and Corwin (2020); Arceri and Corwin (2019).

Here, we move one step further in the comparison between low-dimensional jamming and mean-field jamming by probing not only scaling relations but also prefactors between a handful of properties: pressure $P$, excess contacts $\delta z$, and excess packing fraction above jamming $\Delta\varphi$. We demonstrate the continued success of the mean field in describing low-dimensional systems by quantitatively verifying the mean-field predictions for these prefactors. Thus, the mean-field approximation is overzealous: one need not have vanishing correlations in order to obtain these results. In this spirit we provide a first-principles proof of the relation between pressure and excess packing fraction free of the mean-field assumptions. These results call out for proofs for all of the other universal relations of the jamming transition.

Background – Granular materials undergo a jamming transition at a critical packing fraction $\varphi_{j}$. The number of force bearing contacts between grains jumps abruptly from zero to the minimum number sufficient to support global rigidity and thus global pressure, $Z_{c}$. In a packing of $N$ frictionless, spherical particles in $d$ dimensions, $Z_{c}=Nd+1-d$ Liu and Nagel (2010); Goodrich et al. (2012).

We limit our study to spherical particles interacting through a harmonic contact potential given by

$$U_{ij}=\varepsilon\left(1-\frac{|\mathbf{r}_{ij}|}{\sigma_{ij}}\right)^{2}\Theta\left(1-\frac{|\mathbf{r}_{ij}|}{\sigma_{ij}}\right),$$

(1)

where $\varepsilon$ is the energy scale, $\mathbf{r_{ij}}$ is the contact vector between particles $i$ and $j$, $\sigma_{ij}$ is the sum of the radii of particles $i$ and $j$, and $\Theta$ is the Heaviside step function. Thus, the total energy $U=\frac{1}{2}\sum_{ij}U_{ij}$. From this potential, the forces between particles can be calculated as

$\displaystyle\mathbf{f}_{ij}$

$\displaystyle=\frac{2\varepsilon}{\sigma_{ij}}\left(1-\frac{|\mathbf{r_{ij}}|}{\sigma_{ij}}\right)\Theta\left(1-\frac{|\mathbf{r_{ij}}|}{\sigma_{ij}}\right)\hat{r}_{ij}.$

(2)

We compute a unit and dimension independent pressure using the microscopic formula O’Hern et al. (2003); Allen et al. (1989)

$\displaystyle P$

$\displaystyle\equiv-\frac{\bar{V}_{p}}{\varepsilon}\frac{dU}{dV}=\frac{\bar{V}_{p}}{\varepsilon Vd}\sum_{i,j}\mathbf{f}_{ij}\cdot\mathbf{r}_{ij},$

(3)

where $V$ is the volume of the system and $\bar{V}_{p}$ is the average particle volume.

For soft spheres the packing fraction $\varphi$ can be increased, leading to new contacts and an increased pressure. We thus consider three natural quantities that measure distance from jamming:

• •

  excess packing fraction, $\Delta\varphi=\varphi-\varphi_{j}$

• •

  excess contacts per particle, $\delta z=\left(Z-Z_{c}\right)/N$ where $Z$ is the number of contacts

• •

  pressure $P$

The relationships between these quantities are predicted by mean-field theory as Parisi et al. (2020):

	$\displaystyle P$	$\displaystyle=C_{p\varphi}\Delta\varphi$		(4)
	$\displaystyle\delta z$	$\displaystyle=C_{zp}P^{1/2}$		(5)

with prefactors $C_{p\varphi}$ and $C_{zp}$ which are functions only of spatial dimension O’Hern et al. (2003). These and other scaling relationships have been previously explained by approximate theories Wyart et al. (2005); Wyart (2005); Zaccone and Scossa-Romano (2011); Liarte et al. (2019) and computationally confirmed in low-dimensional jamming O’Hern et al. (2003, 2002); Liu and Nagel (2010); Goodrich et al. (2012). They are summarized concisely by the scaling theory of the jamming transition Goodrich et al. (2016). The scaling exponents in $d\geq 2$ match those in mean field, suggesting that the transition behaves like a critical point with upper critical dimension $d_{u}=2$. Moreover, mean-field theory predictions of these prefactors can be derived as Parisi et al. (2020); Franz et al. (2017):

	$\displaystyle C_{p\varphi}$	$\displaystyle=\frac{1}{d}\hat{C}_{p\varphi}$		(6)
	$\displaystyle C_{zp}$	$\displaystyle=\frac{d}{\sqrt{2^{d}}}\hat{C}_{zp}$		(7)

where $\hat{C}_{p\varphi}$ and $\hat{C}_{zp}$ are finite constants in the $d\rightarrow\infty$ limit, which have not yet been explicitly calculated. Note that these relations are presented in a particular choice of units in the literature. We include details of the conversion to our dimensionless units in the Supplemental Material. A priori, it is not expected that these predictions will apply in low dimensions, in which the mean-field assumption is not warranted. Even above upper critical dimensions, mean-field theories are not generally expected to correctly compute prefactors, or even the purportedly universal amplitude ratios. Beyond scaling exponents, to our knowledge, the critical cluster shape in percolation and related phenomena Aronovitz and Stephen (1987); Privman et al. (1991) and the Binder cumulant in the Ising model Brizin and Zinn-Justin (1985); Parisi and Ruiz-Lorenzo (1996); Blöte and Luijten (1997) are the only quantities which are known to be equal to their mean-field values above the upper critical dimension. Even though these prefactors for jamming scaling relationships have been measured and reported O’Hern et al. (2003); Sartor and Corwin (2020), because they are not expected to be equal to their mean-field values they have not received substantial theoretical attention. An approximate calculation of the related prefactor between the shear modulus and number of excess contacts has been performed in three dimensions Zaccone and Scossa-Romano (2011).

Computational methods – We use pyCudaPacking Charbonneau et al. (2012), a GPU-based simulation engine, to generate energy minimized soft (or penetrable) sphere packings. We do so for number of particles $N=8192-32768$ and dimension $d=2-10$. Our results suggest that $N=8192$ is large enough to avoid finite size effects in $d<9$, which we have verified in $d=8$ by comparing our packing at $N=8192$ with one at $N=16384$, finding no deviation. For $d=9$ and $d=10$ we use system sizes of $16384$ and $32768$, respectively. The particles are monodisperse, except in two dimensions in which we use equal numbers of bidisperse particles with a size ratio of 1:1.4 to prevent crystallization.

The packings are subject to periodic boundary conditions. We minimize the packings using the FIRE minimization algorithm Bitzek et al. (2006) using quad precision floating point numbers in order to achieve resolution on the contact network near the jamming point.

Using the same methods as described in Ref. Charbonneau et al. (2015), we start with randomly distributed initial positions, and apply a search algorithm to create systems approximately logarithmically spaced in $\Delta\varphi$. At each step we use the known power law relationship between energy and $\Delta\varphi$ to calculate an estimate of $\varphi_{j}$. We use this estimate to approximate $\Delta\varphi$ and determine the next value of $\varphi$ in an effort to logarithmically space $\Delta\varphi$ values. We then adjust the packing fraction to this value of $\varphi$ by uniformly scaling particle radii and minimizing the system. We continue this process until the system is nearly critically jammed, i.e. has exactly one state of self stress. We then use the known power law relationship between pressure and $\Delta\varphi$ to fit the dataset and precisely calculate $\varphi_{j}$ (with error less than the smallest value of $\Delta\varphi$) from which we calculate $\Delta\varphi$ at each value of $\varphi$.

Results – Figure 1 shows the measured linear scaling of pressure with packing fraction separately for each dimension. We fit the data to Eq. 4 to find $C_{p\varphi}$, considering only data close to jamming to avoid fitting to high pressure deviations from the scaling power law. The measured values of $C_{p\varphi}$ are shown in the inset to confirm the $1/d$ dimensional scaling predicted by mean-field theory in Eq. 6. A fit to this scaling provides a value of $\hat{C}_{p\varphi}$ of 1.23.

Figure 2 shows the measured square root scaling of excess contacts with pressure separately for each dimension. We fit the data to Eq. 5 to find $C_{zp}$, the values of which are shown in the inset. Beginning around three dimensions, the values of $C_{zp}$ confirm the dimensional scaling predicted by mean-field theory in Eq. 7, and a fit to this scaling provides a value of $\hat{C}_{zp}$ of 0.74.

The values of both $C_{p\varphi}$ and $C_{zp}$ are roughly consistent with values measured in previous studies O’Hern et al. (2003); Sartor and Corwin (2020). It has been recently suggested that the prestress, i.e., the normalized ratio of the first and second derivatives of the potential as defined in Ref. Shimada et al. (2019), is a better candidate to dedimensionalize the relationship between pressure and excess contacts. However, we find a substantially better collapse of our expected form of pressure than with prestress. For more details on prestress, see the attached Supplemental Material.

Discussion – The close agreement of our data with the mean-field predictions in low dimensions suggests that the mean-field assumption is not essential to derive these scaling and prefactor relations. In the spirit of discovering proofs for these relations free of the mean-field assumption, we expand on an earlier calculation of the bulk modulus scaling Wyart (2005) to show that such a calculation can also explain the scaling of $C_{p\varphi}$ with spatial dimension and the precise value of $\hat{C}_{p\varphi}$.

From taking a derivative of Eq. 4, we see immediately that $C_{p\varphi}$ may be expressed in terms of the bulk modulus, $K\equiv V\frac{d^{2}U}{dV^{2}}$, at jamming:

$$C_{p\varphi}=\frac{\bar{V}_{\mathrm{p}}V}{\varphi\varepsilon}\frac{d^{2}U}{dV^{2}}=\frac{V}{N\varepsilon}K.$$

(8)

We note that this approximation slightly overestimates $C_{p\varphi}$: the apparently linear average stress-strain curves of jammed packings are actually the average of many piecewise linear curves with discontinuous drops in stress, thus the average slope is slightly less than the instantaneous slope Fan et al. (2017).

At the unjamming point, the linear response of the system is that of a network of unstretched springs. Thus, at lowest order in pressure the bulk modulus is that of an unstressed spring network, which may be calculated in terms of the “states of self stress,” vectors of possible spring tensions, $s\in\mathbb{R}^{Z}$, which do not produce any net force on a particle Pellegrino (1993); Wyart (2005); Lubensky et al. (2015). Here we explain how to carry out this calculation for a monodisperse system in the unjamming limit; a correction for polydispersity is handled in the Supplemental Material.

We begin by defining the set of “affine bond extensions,” a vector $E\in\mathbb{R}^{Z}$ giving the amount by which each bond vector would increase under a unit volumetric expansion of the system. In linear elasticity, this simply induces an expansion of each length by $1/d$, so,

$$E_{\ell}=\frac{1}{d}r_{\ell},$$

(9)

where we emphasize that $\ell$ indexes the contacts in the system rather than the particles; $r_{\ell}$ is the distance between a particular pair of particles.

In the case that all springs have the same spring constant $k$ (e.g., monodisperse packings), the bulk modulus may be written as the projection of these affine moduli onto the states of self stress Pellegrino (1993); Wyart (2005); Lubensky et al. (2015). At jamming, there is only one state of self stress, and so the bulk modulus may be computed exactly using the projection onto only this one state of self stress Wyart (2005),

	$\displaystyle K$	$\displaystyle=\frac{k}{V}\left(\sum_{\ell=1}^{Z}s_{1,\ell}E_{\ell}\right)^{2}$		(10)
		$\displaystyle=\frac{2N\varepsilon}{dV}\frac{\langle f\rangle^{2}}{\langle f^{2}\rangle}$		(11)

In the near jamming limit, this one special state of self stress exists all the way down to the jamming point and can be expressed in terms of the vector of physical force magnitudes, $f$. For the packing to be in equilibrium, this set of contact forces must produce no net force on every particle, and thus by definition the vector $f$ is always a state of self stress. The projection defined above requires states of self stress to be normalized, and so the state of self stress may be expressed as:

$$s_{1,\ell}=\frac{1}{\sqrt{\sum_{l}f_{l}}}f_{\ell}=\frac{1}{\sqrt{Z\langle f^{2}\rangle}}f_{\ell}.$$

(12)

Furthermore at lowest order in $P$ we have $r=\sigma$, and we assume $Z\approx dN$. Thus, Eq. 10 reduces to

$\displaystyle K$

$\displaystyle=\frac{Nk\sigma^{2}}{dV}\frac{\langle f\rangle^{2}}{\langle f^{2}\rangle}=\frac{2N\varepsilon}{dV}\frac{\langle f\rangle^{2}}{\langle f^{2}\rangle}$

(13)

and thus via Eq. 8

$$C_{p\varphi}=\frac{2}{d}\frac{\langle f\rangle^{2}}{\langle f^{2}\rangle},$$

(14)

for monodisperse spheres. The full calculation in the Supplemental Material shows that in the polydisperse case this becomes

$$C_{p\varphi}=\frac{2}{d}\frac{\langle\sigma f\rangle^{2}}{\langle\sigma^{2}f^{2}\rangle}.$$

(15)

We find that the distribution of contact forces does not depend strongly on dimension, which we demonstrate and discuss in the Supplementary Material, including Refs. Charbonneau et al. (2015); Mueth et al. (1998). We thus predict the scaling of $C_{p\varphi}$ to agree with the asymptotic mean-field scaling. Because this proof does not invoke the mean-field assumption, we expect this scaling to be correct in all dimensions. Moreover, we are able to calculate each value of $C_{p\varphi}$ by measuring the ratio of force distribution moments. These values are calculated as in Eq. 15, and are shown in Fig. 1 to precisely predict the values of $C_{p\varphi}$.

Conclusion – The mean-field theory of jamming predicts both the scaling exponents and the dimensional scaling of their prefactors. While the exponents have been previously verified, we have demonstrated that even some prefactors are well predicted in low dimensions by mean-field theory. Although these prefactors should be considered especially sensitive to finite dimensional corrections, we find the mean field prediction to be exact in low dimensions. Is this a generic phenomenon, or are the quantities we have chosen to study in this work somehow specially unaffected by finite dimensional correlations? Experience with critical phenomena suggests that although certain ratios of these prefactors (i.e. amplitude ratios) may be universal, the prefactors themselves should be both nonuniversal and challenging to compute, which has led to them being neglected. Our results demonstrate however that these prefactors may be computed exactly. These results call out for other theories of jamming and the glass transition which reproduce the mean-field results without such assumptions, or perhaps for a deeper understanding of why certain mean-field computations may be exact in finite dimensions. Additionally, our results suggest that in traditional critical phenomena mean-field theory may compute more for $d\geq d_{u}$ than has been previously appreciated.

Acknowledgments – We thank Francesco Zamponi for generous assistance in interpreting the mean-field results. We also thank Andrea Liu, Jim Sethna, Cam Dennis, and Aileen Carroll-Godfrey for valuable discussion and feedback. This work benefited from access to the University of Oregon high performance computer, Talapas. This work was supported by National Science Foundation (NSF) Career Grant No. DMR-1255370 and the Simons Foundation No. 454939 (J.D.S. and E.I.C.) and by an NSERC PGS-D fellowship and Simons Foundation No. 454945 to Andrea J. Liu. (S.A.R.).