Estimating angoricity and granular equations of state for monodisperse packings of frictional hard spheres in a wide range of densities

In this paper, we study the granular equation of state (EOS) for computer-generated three-dimensional mechanically stable packings of frictional monodisperse particles over a wide range of densities (packing fractions), $\phi =0.56$–0.72. As a statistical physics framework, we utilize the statistical ensemble for granular matter, specifically the “angoricity” ensemble, where the compressional component ${\mathrm{\Sigma }}_p$ of the force-moment tensor serves as granular energy and angoricity $A_p$ is the corresponding granular “temperature.” We demonstrate that the systems under study conform well to this statistical description, and the simple equation of state ${\mathrm{\Sigma }}_p=2.8N{A}_p$ holds very well, where $N$ is the number of particles. We show that granular temperature exhibits a rapid drop around the random-close packing (RCP) limit $\phi \approx 0.64$–0.65, and, hence, one can say that granular packings “freeze” at the RCP limit. We repeat these calculation for shear angoricity $A_{\text{sh}}$ and shear component ${\mathrm{\Sigma }}_{\text{sh}}$ of the force-moment tensor and obtain a similar EOS, ${\mathrm{\Sigma }}_{\text{sh}}=0.85N{A}_{\text{sh}}$. Additionally, we measure the so-called keramicity, an inverse temperature variable corresponding to the determinant of the force-moment tensor, while pressure angoricity corresponds to its trace. We show that inverse keramicity ${\kappa }^{-1}$ and angoricity $A_p$ conform to an EOS $\frac{1}{A_p}\frac{{\mathrm{\Sigma }}_p}{N}+0.11\kappa {\left(\frac{{\mathrm{\Sigma }}_p}{N}\right)}^3=1.2$, whose form is predicted by mean-field theory. Finally, we demonstrate that the alternative statistical ensemble where Voronoi volumes serve as granular energy (and so-called compactivity serves as temperature) does not describe the systems under study well.

Figure 1

Granular energies in the pressure angoricity ensemble per particle ${\varepsilon }_p={\mathrm{\Sigma }}_p/N$ vs packing density $\phi$. Compression components of the force-moment tensor (granular energies) per particle are relatively high at low densities and exhibit a drop at $\phi =0.64$–0.65, around the RCP limit $\phi =0.64$. It happens because $\phi =0.64$ is the lowest density where mechanically stable frictionless packings can exist (and at $\phi =0.65$, crystalline inclusions start to appear). So at $\phi \ge 0.64$, no friction is needed to keep packings stable, and, hence, no normal forces are needed. It means that the force-moment tensor trace can also rapidly drop at $\phi =0.64$–0.65. Error bars indicate $99.7\%$ (3-sigma) confidence intervals, calculated from standard deviation of ${\mathrm{\Sigma }}_p$ for subpackings of size $N_{\text{grp}}=8$.

Figure 2

Local density and local granular energy clouds in the pressure angoricity ensemble. Each point represents a nonoverlapping subpacking of $N_{\text{grp}}=8$ particles along with its local density ${\phi }_{\text{local}}$ and local granular energy per particle ${\varepsilon }_p={\mathrm{\Sigma }}_p/N_{\text{grp}}$ in a given subpacking. It is a visualization of joint probability density functions $f({\phi }_{\text{local}},{\varepsilon }_p)$. There seems to be no strong dependence of a conditional distribution $f\left({\varepsilon }_p\right|{\phi }_{\text{local}})$ on ${\phi }_{\text{local}}$. Red crosses represent global values of each packing as a whole.

Figure 3

Probability density functions (PDFs) for local energies ${\varepsilon }_p={\mathrm{\Sigma }}_p/N_{\text{grp}}$ of subpackings, depending on how many ($M_{\text{grp}}$) subpackings with suitable ${\phi }_{\text{local}}$ (close to packing $\phi$) we select for calculation. Ideally, we need to calculate $f\left({\varepsilon }_p\right|\phi )$, i.e., select only a small number of subpackings whose local densities are close to the global packing density, ${\phi }_{\text{local}}\approx \phi$. This makes histograms very noisy and calculations unreliable. Since $f\left({\varepsilon }_p\right|{\phi }_{\text{local}})\approx f\left({\varepsilon }_p\right)$ (cf. also Fig. 2), we calculate the PDFs by taking more and more subpackings with densities ${\phi }_{\text{local}}$ closest to $\phi$. The plots demonstrate that it is safe to take all the subpackings ($M_{\text{grp}}\approx 1200$) in a given packing (i.e., $f\left({\varepsilon }_p\right)$) to estimate $f\left({\varepsilon }_p\right|\phi )$, which is consistent with Fig. 2.

Figure 4

Log energy differences for the overlapping histograms method [Eq. (13)] in the pressure angoricity ensemble. Here we display the ratio of logarithms of probability density functions for energies per particle from Fig. 3, $f_1\left({\varepsilon }_p\right)/f_2\left({\varepsilon }_p\right)$, along with their fits. We omit from the fitting procedure bins with $<10$ subpackings, and the red dashed lines stop at the edges where at least one of the two overlapping histograms has $<10$ subpackings in a bin. The differences of logarithms look like straight lines indeed, so the method of overlapping histograms is applicable.

Figure 5

Inverse pressure angoricities up to an unknown constant $\frac{1}{A_p}-C_{A,p}$ vs packing density $\phi$. Error bars represent 3-sigma ($99.7\%$) confidence intervals, where standard deviations are taken from 10 estimates of inverse pressure angoricities from histogram overlaps with 10 packings of higher densities.

Figure 6

Pressure angoricities $A_p$ vs packing densities $\phi$. Pressure angoricities exhibit a drop at $\phi \approx 0.64$–0.65, similarly to granular energies ${\mathrm{\Sigma }}_p$ in Fig. 1. Error bars correspond to $95\%$ confidence intervals estimated from the delta method, where we incorporate uncertainty in fits as well as uncertainties from the estimates from Fig. 5.

Figure 7

Granular energy (pressure and shear components of the force-moment tensor) vs different versions of granular temperature (pressure angoricity, shear angoricity, and inverse keramicity). (a) Pressure components of the force-moment tensor per particle $\langle {\varepsilon }_p\rangle =\langle {\mathrm{\Sigma }}_p\rangle /N_{\text{grp}}$ vs angoricities $A_p$ and shear components of the force-moment tensor $\langle {\mathrm{\Sigma }}_{\text{sh}}\rangle /N_{\text{grp}}$ vs shear angoricities $A_{\text{sh}}$ along with respective linear fits $\langle {\mathrm{\Sigma }}_p\rangle =N_{\text{grp}}{k}_{A,p}{A}_p$ and $\langle {\mathrm{\Sigma }}_{\text{sh}}\rangle =N_{\text{grp}}{k}_{A,\mathrm{sh}}{A}_{\text{sh}}$. Here $k_{A,p}=2.81\pm 0.08$ and $k_{A,\mathrm{sh}}=0.85\pm 0.01$ ($95\%$ confidence intervals). (b) Pressure components of the force-moment tensor per particle $\langle {\varepsilon }_p\rangle =\langle {\mathrm{\Sigma }}_p\rangle /N_{\text{grp}}$ vs angoricities $A_p$ when determined from full fits [Eq. (20)] to histogram log-differences, which also account for keramicities $\kappa$. The red line is from the EOS fit in the form of Eq. (22), which can be interpreted as dependence $\langle {\varepsilon }_p\rangle =f(A_p,{\kappa }^{-1})$. The full depiction of the data from simulations and the EOS fit in 3D can be found in Fig. 13 in Appendix pp5. The present panel is the projection of the 3D data onto the plane $A_p,\langle {\varepsilon }_p\rangle$.

Figure 8

Compactivity $\chi$ fits. (a) Average volume of Voronoi cells, $V_{\text{cell}}=\langle V_{\text{grp}}\rangle /N_{\text{grp}}$, where $V_{\text{grp}}$ is the volume of a subpacking of $N_{\text{grp}}=8$ particles. Apart from a small noise due to averaging, it follows an analytical curve $V_{\text{cell}}=\frac{\pi }{6}{D}^3\frac{1}{\phi }$, where the sphere diameter $D=1$. (b) Inverse compactivity $1/\chi$ up to an unknown constant $C_{\chi }$ determined by the method of overlapping histograms. (c) EOS $V_{\text{cell}}$ vs $\chi$ after fitting it with the form $V_{\text{cell}}=a\chi +b/\chi$ (the unknown constant $C_{\chi }$ is fitted implicitly). The red dashed line is the fit. The fit was performed for $\phi \le 0.65$. (d) $\chi$ vs $\phi$. Compare with Fig. 4 from Ref. [17] for 2D. (e) $V_{\text{cell}}\times (1/\chi -C_{\chi })$ vs $V_{\text{cell}}$. If the EOS had the form $V_{\text{cell}}=a\chi$, this plot would be linear, $a-V_{\text{cell}}{C}_{\chi }$. It is not, which means that the EOS is nonlinear. (f) $V_{\text{cell}}/\chi$ vs $\phi$. If the EOS had the form $V_{\text{cell}}=a\chi$, this plot would be a constant, but it is not. Compare with Fig. 2(c) from Ref. [10].

Figure 9

Goodness-of-fit data for log-differences linear fits for compression angoricity (cf. Fig. 4). (a) Log-likelihoods of fits to log energy differences, Eq. (B1). Red dashed line depicts goodness of fit for linear fits and the purple solid line is constructed for third-order polynomial fits for log-differences. The panel shows that third-order polynomials are sufficient at $\phi >0.65$ and for $\phi \le 0.65$ linear fits are as good as third-order fits. At the same time, for pressure angoricity, points at $\phi >0.65$ comply with the same EOS as points $\phi \le 0.65$, so we include them in the final plot Fig. 7 for $A_p$. (b) $p$ values for ${\chi }^2$ statistics calculated for linear fits as a measure of quality of these fits. The panel shows that for none of the densities we can for sure reject the hypothesis that a fit is linear ($p$ values are never below 0.05) but confirms that fit quality deteriorates for $\phi >0.65$.

Figure 10

Details of shear angoricity calculation. Panels (a)–(d) are equivalent to Figs. 1, 3, 4, 5, respectively. Panels (e) and (f) are equivalent to Figs. 9 and 9. (a) Probability distributions of shear components of the force-moment tensor per particle ${\varepsilon }_{\text{sh}}={\mathrm{\Sigma }}_{\text{sh}}/N_{\text{grp}}$. Histograms are calculated for $M_{\text{grp}}\sim 1200$ subpackings (i.e., over $M_{\text{grp}}\sim 1200$ values of ${\varepsilon }_{\text{sh}}$). Red dashed line (circles): $\phi =0.605$; green solid line (diamonds): $\phi =0.641$; and purple solid line (squares): $\phi =0.669$. (b) Log energy differences for compactivity calculations from the overlapping histograms method, Eq. (14). Similarly to the main text, we display here the ratio of logarithms of probability density functions for energies per particle, $f_1\left({\varepsilon }_{\text{sh}}\right)/f_2\left({\varepsilon }_{\text{sh}}\right)$, along with their linear fits. Left blue line: $\phi =0.6430/0.6431$. Right blue line: $\phi =0.6047/0.6082$. (c) Granular energies per particle ${\mathrm{\Sigma }}_{\text{sh}}/N_{\text{grp}}$ vs packing density $\phi$. Error bars have the same meaning as in Fig. 1 (3-sigma confidence intervals). The plot looks qualitatively similar to Fig. 1. (d) Inverse shear angoricities up to an unknown constant $\frac{1}{A_{\text{sh}}}-C_{A,\mathrm{sh}}$ vs packing density $\phi$. Error bars represent 3-sigma ($99.7\%$) confidence intervals, as in Fig. 5. The plot looks qualitatively similar to Fig. 5. (e) Log-likelihoods of fits to log energy differences, Eq. (B1). Red dashed line depicts goodness of fit for linear fits and the purple solid line is constructed for third-order polynomial fits. For $\phi \le 0.65$, linear fits are as good as third-order fits. (f) $p$ values for ${\chi }^2$ statistics calculated for linear fits as a measure of quality of these fits. The panel shows that for none of the densities we can for sure reject the hypothesis that a fit is linear ($p$ values are never below 0.05).

Figure 11

Details of compactivity calculation. Panels (a)–(d) are equivalent to Figs. 3, 4, 9, and 9, respectively. (a) Probability distributions of Voronoi volumes per particle $V_{\text{cell}}=V_{\text{grp}}/N_{\text{grp}}$. Histograms are calculated for $M_{\text{grp}}\sim 1200$ subpackings (i.e., over $M_{\text{grp}}\sim 1200$ values of $V_{\text{grp}}$). Red dashed line (circles): $\phi =0.605$; green solid line (diamonds): $\phi =0.643$; and purple solid line (squares): $\phi =0.669$. (b) Log energy differences for compactivity calculations from the overlapping histograms method, Eq. (14). For $\phi \lesssim 0.65$, the differences of logarithms look like straight lines, so the method of overlapping histograms is applicable. Left blue line: $\phi =0.6685/0.6937$. Middle blue line: $\phi =0.6430/0.6431$. Right blue line: $\phi =0.6047/0.6082$. (e) Log-likelihoods of fits to log energy differences, Eq. (B1). Red dashed line depicts goodness of fit for linear fits and the purple solid line is constructed for third-order polynomial fits for log-differences. For $\phi \le 0.65$, linear fits are as good as third-order fits. Contrary to angoricity plots, even third-order fit quality seems to deteriorate for $\phi >0.65$. (f) $p$ values for ${\chi }^2$ statistics calculated for linear fits as a measure of quality of these fits. The panel shows that for none of the densities we can for sure reject the hypothesis that a fit is linear ($p$ values are never below 0.05).

Figure 12

Investigation of the applicability of using ${\mathrm{\Sigma }}_p^3={\left(\frac{1}{3}\mathrm{Tr}\widehat{\mathrm{\Sigma }}\right)}^3$ instead of $det\widehat{\mathrm{\Sigma }}$, i.e., of fitting histogram differences as one-dimensional curves with Eq. (20) instead of the full fit with a plane $B+\mathrm{\Delta }{\alpha }_p{\mathrm{\Sigma }}_p+\mathrm{\Delta }\kappa det\widehat{\mathrm{\Sigma }}$. For purely isotropic systems, we expect $det\widehat{\mathrm{\Sigma }}=det{\mathrm{\Sigma }}_{p}I={\mathrm{\Sigma }}_p^3$, while for any systems by definition [Eq. (3)] $\frac{1}{3}\mathrm{Tr}\widehat{\mathrm{\Sigma }}={\mathrm{\Sigma }}_p$. The simplified fit with Eq. (20) is applicable only when $det\widehat{\mathrm{\Sigma }}$ is narrowly distributed around ${\mathrm{\Sigma }}_p^3$. The plots demonstrate the applicability of this assumption for three distinct packings with different densities. The panels present heat maps (2D histograms) of subpackings vs ${\mathrm{\Sigma }}_p^3$ and $det\widehat{\mathrm{\Sigma }}$. We take subpackings of eight particles, as used in the main text ($\sim 1200$ subpackings per packing). Colors in each square cell represent the number of subpackings with values of ${\mathrm{\Sigma }}_p^3$ and $det\widehat{\mathrm{\Sigma }}$ inside narrow ranges corresponding to a cell. The colorbar maps colors to a number of particles in a cell. Black dashed lines represent mean values of $det\widehat{\mathrm{\Sigma }}$ vs ${\mathrm{\Sigma }}_p^3$. They are very close to linear lines. Blue dots follow the line $det\widehat{\mathrm{\Sigma }}={\left(\frac{1}{3}\mathrm{Tr}\widehat{\mathrm{\Sigma }}\right)}^3$.

Figure 13

Fitting histogram log-differences with higher-order fits (estimating keramicity), cf. Ref. [12] and Eq. (19). (a) Inverse pressure angoricities up to a constant when simultaneously estimating keramicities. This panel is equivalent to Fig. 5 but estimated from Eq. (20). The panel is qualitatively similar to Fig. 5. (b) Keramicity, or inverse “Beltrami volume temperature” [cf. Eq. (20)] up to an unknown constant. The error bars in panels (a) and (b) represent 3-sigma ($97\%$) confidence intervals. (c) Complete angoricity-keramicity EOS [Eq. (22)] estimated from inverse angoricity and keramicity values from panels (a) and (b). Blue dots depict triples $(A_p,{\kappa }^{-1},\langle {\varepsilon }_p\rangle )$, where ${\varepsilon }_p={\mathrm{\Sigma }}_p/N_{\text{grp}}$, determined from simulations for each granular packing. The surface represents the EOS from fitting Eq. (22), which can be interpreted as a 3D plot $\langle {\varepsilon }_p\rangle =f(A_p,{\kappa }^{-1})$. Figure 7 is a projection of this panel onto the axes $A_p,\langle {\varepsilon }_p\rangle$. (d) Inverse pressure angoricities up to a constant from a linear fit (13), keramicity not included, vs granular energy per particle ${\varepsilon }_p={\mathrm{\Sigma }}_p/N_{\text{grp}}$. Red circles are from packings with $\phi \le 0.644$; blue circles are points with $\phi >0.644$. The plot is a hyperbola and can indeed be fitted by Eq. (16). It shows that one can fit the data at all densities with the same Eq. (16), which we do in the main text, Fig. 7. (e) Same as panel (d) but with inverse angoricities from panel (a). (f) Keramicities up to a constant from panel (b) vs granular energy per particle ${\varepsilon }_p={\mathrm{\Sigma }}_p/N_{\text{grp}}$. Red circles are from packings with $\phi \le 0.644$; blue circles are points with $\phi >0.644$. Points from packings with $\phi \le 0.644$ seem to lie on the same curve, meaning they can be fitted with Eq. (22). The plot shows that packings with $\phi >0.644$ do not comply with a EOS, presumably because keramicity is more sensitive to packings having a random structure, so we exclude them from fitting keramicity EOS [cf. panel (c) and Fig. 7].

Figure 14

Analysis of finite-size effects on the granular EOS for pressure angoricity. We repeated the measurements from the main text to estimate $k_{A,p}$ from Eq. (16) for different group sizes of subpackings, $N_{\text{grp}}=7$, 8, 12, 16, and 32. All the results for different group sizes look qualitatively similar to the results presented in the main text. We fitted this curve in the form $k_{A,p}=f(1/N_{\text{grp}}^2)$, where $f$ is the second-order polynomial. The estimate of $k_{A,p}$ in the thermodynamic limit $N_{\text{grp}}\to \infty$ is $k_{A,p}=2.60\pm 0.27$. The error bars represent $95\%$ confidence intervals.