Jamming in finite systems: Stability, anisotropy, fluctuations, and scaling

Athermal packings of soft repulsive spheres exhibit a sharp jamming transition in the thermodynamic limit. Upon further compression, various structural and mechanical properties display clean power-law behavior over many decades in pressure. As with any phase transition, the rounding of such behavior in finite systems close to the transition plays an important role in understanding the nature of the transition itself. The situation for jamming is surprisingly rich: the assumption that jammed packings are isotropic is only strictly true in the large-size limit, and finite-size has a profound effect on the very meaning of jamming. Here, we provide a comprehensive numerical study of finite-size effects in sphere packings above the jamming transition, focusing on stability as well as the scaling of the contact number and the elastic response.

Figure 1

Schematic energy landscape where $|r\rangle$ denotes the particle degrees of freedom and $|\Delta L\rangle$ all possible shear deformations of the box. For packings in the ${\mathcal{E}}_{\mathrm{comp}}$ ensemble, $|\Delta L\rangle$ is fixed and the system is ${\mathcal{R}}_{\mathrm{comp}}$ stable if it sits at a minimum of $U$ with respect to $|r\rangle$, i.e., the open circle. ${\mathcal{R}}_{\mathrm{all}}$ stability is governed by the curvature of $U$ along the global shear degrees of freedom. Thus, ${\mathcal{R}}_{\mathrm{comp}}$ stable states can be ${\mathcal{R}}_{\mathrm{all}}$ unstable if the curvature of $U$ is negative along any of the $|\Delta L\rangle$ directions (thick dashed curve). Such states can and do occur. If the curvature of $U$ is positive along all global shear directions (thick solid curve), the packing is ${\mathcal{R}}_{\mathrm{all}}$ stable. Such ${\mathcal{R}}_{\mathrm{all}}$ stable packings can have finite shear stresses (nonzero gradient along global shear directions). Finally, packings that are at a local minimum of $U$ with respect to the $|r\rangle$ and $|\Delta L\rangle$ directions (filled circle) have zero residual shear stress in addition to being ${\mathcal{R}}_{\mathrm{all}}$ stable and thus satisfy the ${\mathcal{R}}_{\mathrm{all}}^{+}$ requirement.

Figure 2

Schematic depiction of the hierarchical relation among the three ensembles ${\mathcal{E}}_{\mathrm{comp}}$, ${\mathcal{E}}_{\mathrm{all}}$, and ${\mathcal{E}}_{\mathrm{all}}^{+}$, for small $p{N}^2$ (left) and for large $p{N}^2$ (right). Whereas for $p{N}^2\to 0$, ${\mathcal{E}}_{\mathrm{all}}$ becomes vanishingly small in comparison to ${\mathcal{E}}_{\mathrm{comp}}$ [16], in the thermodynamic limit, ${\mathcal{E}}_{\mathrm{comp}}$ and ${\mathcal{E}}_{\mathrm{all}}$ become virtually indistinguishable for all finite $p$. Moreover, for large systems the ratio of residual shear stress to pressure vanishes, so the properties of ${\mathcal{E}}_{\mathrm{all}}$ and ${\mathcal{E}}_{\mathrm{all}}^{+}$ converge.

Figure 3

(a) Example of the sinusoidal function $G\left(\theta \right)$ for a two-dimensional system of $N=256$ particles at a pressure of $p={10}^{-3}$. The black open circles show direct numerical calculations using various $\epsilon \left(\theta \right)$, while the solid line shows the prediction using the elastic modulus tensor at $\theta =0$ (see Appendix pp2). The system is stable to shear when $G\left(\theta \right)>0$ and unstable when $G\left(\theta \right)\le 0$. The data agree with the prediction. The horizontal dashed line shows the average $G_{\mathrm{DC}}$; $G_{\mathrm{AC}}$ is obtained from the amplitude of the sinusoidal curve.

Figure 4

The total number of contacts $N_c$ above the isostatic number $N_c^{\mathrm{iso}}$ as a function of pressure for systems of $N=256$ particles. The solid horizontal line is at $N_c-N_c^{\mathrm{iso}}=1$ and the dashed horizontal line is at $N_c-N_c^{\mathrm{iso}}=3$.

Figure 5

The excess contact number ($Z-Z_{\mathrm{iso}}^N$, top row), shear modulus ($G_{\mathrm{DC}}$, middle row), and bulk modulus ($B$, bottom row) as a function of pressure for system sizes ranging from $N=32$ (blue squares) to $N=4096$ (red crosses). The four columns show, from left to right, the $2d$ ${\mathcal{E}}_{\mathrm{comp}}$ ensemble, the $3d$ ${\mathcal{E}}_{\mathrm{comp}}$ ensemble, the $2d$ ${\mathcal{E}}_{\mathrm{all}}^{+}$ ensemble, and the $2d$ ${\mathcal{E}}_{\mathrm{all}}$ ensemble. The solid black lines all have slope $1/2$.

Figure 6

Finite-size scaling collapse of the excess contact number and average elastic moduli. Top row: $Z$ minus the theoretical minimum $Z_{\min }^N$ [see Eq. (19)]. $(Z-Z_{\min }^N)N$ collapses as a function of $p{N}^2$ for the $2d$ ${\mathcal{E}}_{\mathrm{comp}}$ (left), $3d$ ${\mathcal{E}}_{\mathrm{comp}}$ (middle), and $2d$ ${\mathcal{E}}_{\mathrm{all}}^{+}$ (right) ensembles. Note that at low $p{N}^2$, most of the ${\mathcal{E}}_{\mathrm{comp}}$ packings are ${\mathcal{R}}_{\mathrm{all}}$ unstable and our filtered, $2d$ ${\mathcal{E}}_{\mathrm{all}}$ ensemble does not have many states at low $p{N}^2$. This is why data are not shown for this ensemble. At large $p{N}^2$, $Z-Z_{\min }^N\sim Z-Z_{\mathrm{iso}}^N\sim p^{1/2}$ [Eq. (20)], while $Z-Z_{\min }^N\sim pN$ at low $p{N}^2$ [Eq. (23)]. The crossover between these scalings occurs when the total number of extra contacts is of order 10. Middle row: $G_{\mathrm{DC}}$ minus the measured $p\to 0$ plateau. $(G_{\mathrm{DC}}-G_{DC,0})N$ collapses as a function of $p{N}^2$ and has the same crossover behavior as $(Z-Z_{\min }^N)N$. The insets show that $G_{DC,0}$ is proportional to $N^{-1}$. Bottom row: $B$ minus the measured $p\to 0$ plateau. Note that the plateau $B_0$ of the bulk modulus is much larger than for the shear modulus. Therefore, uncertainties in $B$ lead to the large error bars in $(B-B_0)N$ at low $p{N}^2$. The insets show that $B_0$ is roughly constant in $N$, as expected for particles with harmonic interactions. It is not clear from the data whether there is an additional $N^{-1}$ contribution to the plateau [i.e., $B_0\left(N\right)=B_0\left(\infty \right)+a{N}^{-1}$]. The colors and symbols are the same as in Fig. 5.

Figure 7

(a) $(Z-Z_{\min }^N)p^{-1/2}$ as a function of $p{N}^2$ on a linear scale for the two- and three-dimensional ${\mathcal{E}}_{\mathrm{comp}}$ ensembles. The $3d$ data show good collapse but there is a system size dependence in $2d$ that was not as clear in Fig. 6. The two-dimensional ${\mathcal{E}}_{\mathrm{all}}^{+}$ data (not shown) is indistinguishable from the $2d$ ${\mathcal{E}}_{\mathrm{comp}}$ data. (b) $(Z-Z_{\min }^N)p^{-1/2}{\left({log}_{10}N\right)}^{y/2}$ as a function of $p{N}^2{\left({log}_{10}N\right)}^{-y}$, with $y=0.7$. The $2d$ ${\mathcal{E}}_{\mathrm{comp}}$ and $2d$ ${\mathcal{E}}_{\mathrm{all}}^{+}$ ensembles are both shown and collapse perfectly onto each other. This shows that the system size dependence in (a) can be accounted for by introducing a logarithmic correction of the form of Eq. (25). The colors and symbols are the same as in Fig. 5.

Figure 8

The average “$\mathrm{AC}$” quantities, which are defined in Appendix pp2 and discussed in the text. The colors and symbols are the same as in Fig. 5.

Figure 9

Scaling collapse of $U_{\mathrm{AC}}$ and $D_{\mathrm{AC}}$ for the $2d$ ${\mathcal{E}}_{\mathrm{all}}^{+}$ ensemble. The scaling of these two quantities is the only unexpected difference between the three ensembles that we have found. The colors and symbols are the same as in Fig. 5.

Figure 10

(a) The probability distribution $P(s/〈s〉)$ of the residual shear stress divided by the ensemble average. $P(s/〈s〉)$ collapses onto a single curve and is independent of system size and pressure. Inset: $〈s〉N^{1/2}$ as a function of pressure. (b) $U_{\mathrm{AC}}$ as a function of $p{N}^2$ for $N=256$ systems from the $2d$ ${\mathcal{E}}_{\mathrm{comp}}$ ensemble. Systems are binned according to the residual shear stress before averaging. For high $s$ (green diamonds), $U_{\mathrm{AC}}$ is roughly constant but for low $s$ (blue squares), $U_{\mathrm{AC}}$ is similar to the ${\mathcal{E}}_{\mathrm{all}}^{+}$ data in Fig. 9. (c) $D_{\mathrm{AC}}$ displays the same behavior as $U_{\mathrm{AC}}$. [(d) and (e)] Scatter plot of the lowest pressure values of $U_{\mathrm{AC}}$ and $D_{\mathrm{AC}}$, both of which show a remarkable linear dependence on the shear stress. The colors and symbols in (a), (d), and (e) are the same as in Fig. 5.

Figure 11

The relative fluctuations in $Z$, $G_{\mathrm{DC}}$, and $B$. The colors and symbols are the same as in Fig. 5.

Figure 12

The distribution of $G_{\mathrm{DC}}$ for the different ensembles. The shape of the distribution function differs for low $p{N}^2$ (top row), medium $p{N}^2$ (middle row), and high $p{N}^2$ (bottom row) but collapses for systems at similar $p{N}^2$. The precise values of $p{N}^2$ correspond to the vertical dashed lines in Fig. 11 (${10}^{-1}$, ${10}^0$, and ${10}^1$ for the $2d$ ensembles; ${10}^{-1.5}$, ${10}^{-0.5}$ and ${10}^{0.5}$ for the $3d$ ensemble).