Critical scaling for yield is independent from distance to isostaticity

TL;DR

This study shows that the critical behavior of yielding in soft particle packings is independent of how close the system is to the isostatic point, revealing a distinct nonequilibrium critical transition.

Contribution

It demonstrates that the critical scaling for yielding is universal across different pressures, independent of the isostaticity proximity, contrasting with previous isostatic critical transition studies.

Findings

Critical scaling functions are nearly independent of pressure.

Diverging length scale follows a9 a9 |a9-a9_c|^{- u}.

Yielding is a distinct nonequilibrium critical transition.

Abstract

Using discrete element simulations, we demonstrate that critical behavior for yielding in soft disk and sphere packings is independent of distance to isostaticity over a wide range of dimensionless pressures. Jammed states are explored via quasistatic shear at fixed pressure, and the statistics of the dimensionless shear stress $\mu$ of these states obey a scaling description with diverging length scale $\xi \propto |\mu-\mu_c|^{-\nu}$. The critical scaling functions and values of the scaling exponents are nearly independent of distance to isostaticity despite the large range of pressures studied. Our results demonstrate that yielding of jammed systems represents a distinct nonequilibrium critical transition from the isostatic critical transition which has been demonstrated by previous studies. Our results may also be useful in deriving nonlocal rheological descriptions of granular materials, foams, emulsions, and other soft particulate materials.

Equations

$${\sigma }_{\alpha \lambda }=\frac{1}{V}\sum _{i\ne j}{r}_{\alpha }^{ij}{F}_{\lambda }^{ij}.$$

$$F$$

$$⟨{\gamma }_{\mathrm{slip}}⟩$$

$$⟨{\gamma }_{\mathrm{ms}}⟩$$

Full text

Granular materials, dense suspensions, foams, and emulsions can form amorphous jammed states Jaeger et al. (1996); O’Hern et al. (2003); Donev et al. (2004); Van Hecke (2009); Liu and Nagel (2010). Jammed systems can yield when subjected to a sufficient shear stress $\tau$ (this is sometimes called “unjamming by shear”). When $\mu\equiv\tau/p$, where $p$ is the system pressure, exceeds a critical value $\mu_{c}$, jammed states become inaccessible and flow persists indefinitely Drucker and Prager (1952); Gardiner et al. (1998); Toiya et al. (2004); da Cruz et al. (2005); Xu and O’Hern (2006); Jop et al. (2006); Peyneau and Roux (2008a, b). Predicting the mechanical response of jammed states can be difficult since it can involve plastic rearrangement events that cooperate over large distances. For example, rheological models of these materials that include nonlocal cooperative effects can successfully reproduce steady-state flows from experiments and particle-based simulations Goyon et al. (2008); Masselon et al. (2008); Bocquet et al. (2009); Jop et al. (2012); Kamrin and Koval (2012); Henann and Kamrin (2013); Bouzid et al. (2013, 2015).

Prior studies on soft sphere packings, which are commonly used to model these materials, have framed long-range cooperative behavior in terms of a nonequilibrium critical transition that occurs at the isostatic point, also called “point J” O’Hern et al. (2003); Drocco et al. (2005); Ellenbroek et al. (2006); Olsson and Teitel (2007); Nordstrom et al. (2010); Tighe et al. (2010); Olsson and Teitel (2011); Paredes et al. (2013). Isostaticity refers to the number of contacts per particle $Z$ being equal to the number required to constrain all degrees of freedom in the system, $Z=Z_{\rm iso}$. This occurs at a given volume fraction $\phi=\phi_{c}$ in the large-system limit. At isostaticity, $p=0$, but further compression leads to increasing $p$. A cooperative length scale $\xi_{J}\propto|\phi-\phi_{J}|^{-\nu_{J}}$ diverges at the isostatic point, which then controls the mechanical response and leads to Widom-like scaling relations Goodrich et al. (2016) that relate $p$, $(\phi-\phi_{J})$, $(Z-Z_{\rm iso})$, and other quantities. $\xi_{J}$ is large near isostaticity (i.e., small $p$), characterized by an excess of spatially extended, low-energy modes of the system Silbert et al. (2005); Wyart et al. (2005). For increasing $p$, $\xi_{J}$ decreases, leading to more localized modes as well as smaller and more localized particle rearrangements.

In contrast, nonlocal rheological descriptions of jammed materials Bocquet et al. (2009); Kamrin and Koval (2012); Henann and Kamrin (2013); Bouzid et al. (2013, 2015); Clark et al. (2018) often include a diverging cooperative length scale that depends not on packing fraction but on distance to a critical shear stress, i.e., $\xi\propto|\mu-\mu_{c}|^{-\nu}$. These rheological models, including our previous paper Clark et al. (2018), describe materials that are near $\phi=\phi_{J}$, so it is not known how the cooperative length scale underlying these models relates to the isostatic critical point. Here we show using numerical simulations that yielding in soft sphere packings is a distinct nonequilibrium critical transition and that it is independent from distance to isostaticity. We quasistatically shear systems of repulsive, bidisperse disks and spheres, holding dimensionless pressure $\tilde{p}$ fixed and measuring $\mu$, which increases during an initial shear regime and then plateaus as stress is released in intermittent slips. The statistics of $\mu$ obey a scaling description with a diverging length scale $\xi\propto|\mu-\mu_{c}|^{-\nu}$, where $\nu_{\rm ms}\approx 1.8$ during initial shear buildup (in agreement with Clark et al. (2018)) and $\nu_{\rm slip}\approx 1.1$ in two dimensions (2D) and $\nu_{\rm slip}\approx 0.8$ in three dimensions (3D) during the intermittent slip regime. The scaling functions and the values of $\nu$ are highly insensitive to the distance from isostaticity set by $\tilde{p}$, which we vary over nearly four orders of magnitude. $\mu_{c}(\tilde{p})$ is insensitive to $\tilde{p}$ for $\tilde{p}\leq 10^{-3}$, but decreases logarithmically for higher $\tilde{p}$. The critical scaling functions we show could be used to derive a particle-scale theory for nonlocal rheological models, including transient behavior, which is not captured by current models.

Methods.— We use molecular dynamics simulations to study systems of $N$ bidisperse frictionless disks in 2D and spheres in 3D, with diameter ratio $1.4$ and equal numbers of each size. Systems are prepared at a given pressure $p$ via isotropic compression and then quasistatically sheared. Contacting particles interact via a purely repulsive force $\mathbf{F}_{ij}=K(\delta_{ij}/|\mathbf{r}_{ij}|-1)\mathbf{r}_{ij}$, where $\delta_{ij}$ is the average diameter of particles $i$ and $j$, $\mathbf{r}_{ij}$ is the vector displacement between the centers of particles $i$ and $j$. Stresses are quantified by the Cauchy stress tensor,

[TABLE]

Here, $\alpha$ and $\lambda$ are Cartesian coordinates, $V$ is the system volume, $r^{ij}_{\alpha}$ is the $\alpha$-component of the center-to-center separation vector between particles $i$ and $j$, and $F^{ij}_{\lambda}$ is the $\lambda$-component of the interparticle contact force. The sum over $i$ and $j$ includes all pairs of contacting particles.

Each simulation is prepared by placing particles randomly throughout the domain and then increasing the particle diameter $D$ in small steps until reaching a target $p=(\sigma_{xx}+\sigma_{yy})/2$. Using Lees-Edwards boundary conditions, we impose affine shear strain in small steps $\Delta\gamma=10^{-4}$. At each shear step, the shear-periodic boundary is moved by $\Delta\gamma$ and $y\Delta\gamma$ is added to the $x$-position of each particle. We then use molecular dynamics to relax the potential energy using modified velocity Verlet integration, as well as shrink or swell $D$ to maintain a fixed $p$ within $0.5\%$ of the target value. Before shearing or changing the particle diameter, we damp out kinetic energy via a viscous damping force $-B\mathbf{v}$ to each particle, where $\mathbf{v}$ is the absolute velocity of a given particle and $B$ is the damping coefficient. We set $B=5\sqrt{Mp}$, where $M$ is the mass of a large grain. Our results are independent of $B$ in this regime.

At each strain step, after the system is quenched at the target pressure, we measure the stress tensor elements, as defined in Eq. (1), focusing on $\mu=-\sigma_{xy}/p$, as shown in Fig. 1(c). We measure $\mu$ from $0\leq\gamma\leq 3$ in increments of $\Delta\gamma=10^{-4}$ for a total of 30,001 states per simulation. For each value of $N$ and $\tilde{p}$, we simulate an ensemble of 400 systems.

We quantify distance above isostaticity by $\tilde{p}=p/K$, which gives an estimate of the relative overlap between particles (i.e., $\tilde{p}=0.01$ corresponds to particle-particle overlaps of roughly $0.01D$). Figure 1(a) and (b) show $\tilde{p}=0.2$ and $\tilde{p}=0.001$, respectively, with $N=24$. Overcompression yields excess contacts such that $\Delta Z/Z_{\rm iso}\sim\sqrt{\tilde{p}}$, where $\Delta Z\equiv(Z-Z_{\rm iso})$, $Z$ is the number of contacts per particle, $Z_{\rm iso}\equiv 2\left(d\left(N-N_{r}\right)-d+1\right)/\left(N-N_{r}\right)$ is the isostatic number of contacts, $d$ is the number of spatial dimensions, and $N_{r}$ is the number of rattlers Van Hecke (2009); Papanikolaou et al. (2013); Shen et al. (2014); VanderWerf et al. (2018). Figure 2(a) shows $\Delta Z/Z_{\rm iso}$ plotted versus $\gamma$ for $N=50$ and varying $\tilde{p}$. These curves fluctuate around a fixed value but show no trend, indicating that the shearing does not change $Z$ on average. Figure 2(b) shows that the average value $\langle\Delta Z\rangle/Z_{\rm iso}$ versus $\tilde{p}$ is similar for $N=50$, $100$, and $200$. Thus, the fraction of excess contacts and thus the distance to isostaticity is set by $\tilde{p}$, nearly independent of system size Goodrich et al. (2012) or the presence of shear deformation Favier de Coulomb et al. (2017).

Scaling near yield.— As shown in Fig. 1(c), $\mu$ increases with $\gamma$ and then plateaus as potential energy is released in intermittent slips Miguel et al. (2001); Dalton and Corcoran (2001); Maloney and Lemaître (2006); Daniels and Hayman (2008); Dahmen et al. (2011); Salerno et al. (2012); Barés et al. (2017). This curve represents a series of jammed states that the system passes through while sheared. The fluctuations in $\mu$ decrease with increasing $N$ for all $\tilde{p}$, and we exploit the size scaling in these fluctuations (as in Lin et al. (2014)) to demonstrate and quantify diverging spatial correlations near $\mu_{c}$. Most importantly, we show that this scaling description is nearly independent from the distance to isostaticity.

To accomplish this, we use finite size scaling on three quantities for each $\mu$ and $N$: (1) the cumulative distribution function $F$ of states above a particular value of $\mu$ during the slip avalanche regime, defined as $\gamma>0.5$ (our results are insensitive to this choice); (2) the shear strain $\gamma_{\rm slip}$ between mechanically stable (MS) states with an internal shear stress of at least $\mu$; and (3) the shear strain $\gamma_{\rm ms}$ required to find the first MS state with an internal shear stress of at least $\mu$. Figure 1(c) depicts $\gamma_{\rm ms}$ and $\gamma_{\rm slip}$ for a given $\mu(\gamma)$ curve. Figure 3(a)-(c) shows these quantities plotted as functions of $\mu$ or $N$, where ensemble averages are denoted with angle brackets. The data shown in Fig. 3 represents only a single value of $\tilde{p}=0.05$ in 2D, but it is typical of all $\tilde{p}$ in both 2D and 3D, as we demonstrate below in Figs. 4 and 5. As $N$ is increased, the fluctuations in $\mu$ decrease, and $F$ approaches a step function, as shown in Fig. 3(a). Thus, MS states vanish sharply at some value $\mu=\mu_{c}(\tilde{p})$ in the large system limit. Figure 3(b) shows $\langle\gamma_{\rm slip}\rangle$, where we require at least one $\gamma_{\rm slip}$ measurement per simulation. Our results are insensitive to this choice, unless the number of samples becomes very small. For $\mu<\mu_{c}$, $\langle\gamma_{\rm slip}\rangle$ monotonically decreases with increasing $N$. For $\mu>\mu_{c}$, $\langle\gamma_{\rm slip}\rangle$ first decreases and then increases with increasing $N$. Finally, $\langle\gamma_{\rm slip}\rangle$ is nearly independent of $N$ for small $\mu$ and increases strongly with $N$ for larger $\mu$.

To collapse these curves, we posit a diverging length scale $\xi\propto\left|\mu-\mu_{c}\right|^{-\nu}$. In this case, finite size effects should enter through the quantity $L/\xi$, where $L\equiv N^{1/d}$ with $d$ being the number of spatial dimensions. An equivalent scaling can also be written using $\left(\mu-\mu_{c}\right)L^{1/\nu}$; see Refs. Olsson and Teitel (2011); Clark et al. (2018) for further discussion on similar systems. Figure 3(d)-(f) shows that the data in Fig. 3(a)-(c) collapses according to

[TABLE]

Here, $f_{\rm ms,\pm}$ and $f_{\rm slip,\pm}$ are dual-branch functions, with $+$ and $-$ denoting $\mu$ above or below $\mu_{c}$, respectively. Interestingly, we need two distinct values of $\nu$ for the initial strain, $\nu_{\rm ms}\approx 1.8$, and slip avalanche regime, $\nu_{\rm slip}\approx 1.1$ in 2D or $\nu_{\rm slip}\approx 0.8$ in 3D. The value $\nu_{\rm slip}\approx 1.8$ agrees with our previous result Clark et al. (2018), which was only calculated near to isostaticity; in Fig. 3, it is calculated far from isostaticity. The difference between $\nu_{\rm slip}$ and $\nu_{\rm ms}$ suggests that there are important differences in how MS states are accessed between these two regimes. We obtain the critical parameters $\mu_{c}$ and $\nu_{\rm slip}$ by fitting the collapsed curves to appropriate functional forms using a Levenberg-Marquardt method Olsson and Teitel (2011); Clark et al. (2018). We exclude system sizes with $N<N_{\rm min}$, and vary $N_{\rm min}$ until our fits become insensitive to our choice of $N_{\rm min}$. We also performed the corrections-to-scaling analysis in Vågberg et al. (2011), which yields the same result we find with the scaling forms in Eqs. (2)-(4).

We then perform the same analysis for varying $\tilde{p}$ over the ranges $5\cdot 10^{-5}\leq\tilde{p}\leq 2\cdot 10^{-1}$ in 2D and $2\cdot 10^{-4}\leq\tilde{p}\leq 2\cdot 10^{-1}$ in 3D, spanning from near isostaticity (where $\xi_{J}$ is large) to far from isostaticity (where $\xi_{J}$ is small). The scaling description in Eqs. (2)-(4) and shown in Fig. 3 holds for all values of $\tilde{p}$ in both 2D and 3D. We show data for an additional pressure in 2D, $\tilde{p}=0.0005$, in Figure 4. We also show data in 3D with $\tilde{p}=0.05$ in Fig. 5. In both cases, the scaling functions are almost indistinguishable from those shown in Fig. 3. Figure 6 shows the critical parameters $\mu_{c}$ and $\nu_{\rm slip}$ plotted as a function of $\tilde{p}$. Each data point in Fig. 6 represents a fit of all data (as in Figs. 3 through 5) over many system sizes (typically $32\leq N\leq 1024$) with 400 simulations per system size, so the plateau in Fig. 6 is not a system size effect. We find $\mu_{c}$ to be independent of $\tilde{p}$ for $\tilde{p}\leq 10^{-3}$, and $\mu_{c}$ decreases logarithmically for $\tilde{p}>10^{-3}$, which agrees with Favier de Coulomb et al. (2017). This occurs as excess contacts are added, which changes the structure of the force networks.

However, we observe no similar crossover behavior as distance to isostaticity is varied in any other aspects of the scaling behavior. The critical exponents, as shown in Fig. 6(b), and the scaling functions, as shown in Figures 3, 4, and 5, are highly insensitive to $\tilde{p}$, despite the wide variation in distance to isostaticity. Specifically, we find $\nu_{\rm slip}\approx 1.1\pm 0.1$ in 2D, $\nu_{\rm slip}\approx 0.8\pm 0.03$ in 3D, $\beta_{\rm slip}/\nu_{\rm slip}\approx-1\pm 0.1$ in 2D, and $\beta_{\rm slip}/\nu_{\rm slip}\approx-1.3\pm 0.1$ in 3D. The uncertainty is estimated from the scatter in the data for different $\tilde{p}$, as seen in Fig. 6. For the initial shear regime, we find that $\nu_{\rm ms}\approx 1.8$ is insensitive to $\tilde{p}$. However, $\beta_{\rm ms}/\nu_{\rm ms}$ appears to vary from roughly 0.2 at high $\tilde{p}$ to 0.6 at low $\tilde{p}$. This again points to potentially important differences between how MS states are explored between the slip avalanche and initial shear regimes and may have consequences for size-dependent arrest transitions Kamrin and Henann (2015); Clark et al. (2018); Srivastava et al. (2018).

Discussion.— We have shown here that sheared amorphous soft sphere packings display finite size scaling that is consistent with a diverging length scale $\xi\propto|\mu-\mu_{c}|^{-\nu}$. The value of $\mu_{c}$ varies as $\tilde{p}$ is changed and extra contacts are added, but the forms of the scaling functions (as shown in Figures 3 , 4, and 5) and the values of the critical exponents are nearly independent of distance to isostaticity over nearly four orders of magnitude in $\tilde{p}$. Considering the correlation length for isostaticity $\xi_{J}\propto|\phi-\phi_{c}|^{-\nu_{J}}$, if one assumes that $\nu_{J}$ is order unity Olsson and Teitel (2011) and $\tilde{p}\propto(\phi-\phi_{c})$ for harmonic interactions O’Hern et al. (2003), then varying $\tilde{p}$ over this range represents $\xi_{J}$ varying over a similar range. This represents an enormous variation with respect to the isostatic critical point, implying that the distance to isostaticity does not control the critical behavior we demonstrate here. Our results suggest that yielding in, e.g., emulsions, foams, or granular materials is controlled by an underlying nonequilibrium critical transition that is distinct from isostaticity. We note that $\nu_{\rm slip}\approx 1.1$ for 2D and $\nu_{\rm slip}\approx 0.8$ in 3D are similar to the values $\nu\approx 1.1$ for 2D and $\nu\approx 0.7$ for 3D from Ref. Lin et al. (2014), which presented a scaling description for yielding in amorphous materials Miguel et al. (2001); Maloney and Lemaître (2006); Salerno et al. (2012).

Figure 6(c) shows the Liu-Nagel jamming phase diagram from, e.g., Refs. O’Hern et al. (2003); Heussinger and Barrat (2009); Bi et al. (2011) and many others, but with $\tilde{p}$ on the horizontal axis and $\mu=\tau/p$ on the vertical axis. The solid blue line represents the critical yielding boundary in 2D from Fig. 6(a), and the solid vertical black line represents the isostatic critical transition. Jammed states exist only in the lower right region, above isostaticity and below the critical yielding boundary. Unjammed or fluid-like states can be either hypostatic ($Z<Z_{\rm iso}$ and $p=0$) or hyperstatic ($Z>Z_{\rm iso}$ and $p>0$). Some previous work on critical scaling near isostaticity has studied the onset of yield stress behavior under shear at varied $\phi$ Olsson and Teitel (2007); Heussinger and Barrat (2009); Nordstrom et al. (2010); Olsson and Teitel (2011); Paredes et al. (2013). Such a system is situated at the “triple point” indicated by the red dot at the intersection of the jamming and yielding lines in Fig. 6(c). A complete theory may be able to unify these two critical transitions by a better understanding of the behavior at this point.