Yielding versus jamming: critical scaling of sheared, soft-core disks and spheres

Using discrete element simulations, we demonstrate critical behavior for yielding of assemblies of soft-core repulsive disks (2D) and spheres (3D) over a wide range of dimensionless pressures. Assemblies are isotropically compressed, and we then perform quasi-static simple shear at fixed pressure using shear-periodic boundaries. By examining the fluctuations in the dimensionless shear stress $\Sigma$, we observe finite-size scaling consistent with a diverging length scale $\xi \propto |\Sigma - \Sigma_c|^{-\nu}$. We observe two distinct values of $\nu$: $\nu_{\rm ms} \approx 1.8$ characterizes the initial stress buildup in both 2D and 3D, and $\nu_{\rm slip}$ characterizes slips during steady-state shear, where $\nu_{\rm slip}\approx 1.1$ in 2D and $\nu_{\rm slip}\approx 0.8$ in 3D. The critical stress $\Sigma_c$ is constant for low pressure, $\Sigma_c \approx 0.1$ in 2D and $\Sigma_c \approx 0.11$ in 3D, but decreases for larger pressures. However, the critical behavior, including the values of scaling exponents, is otherwise unchanged over a wide range of pressures, including far from the jamming transition. Our results show that yielding is in fact a distinct phase transition from jamming, which may explain similarities between nonlocal rheological descriptions of granular materials, foams, emulsions, and other soft particulate materials.

Granular materials, dense suspensions, foams, emulsions, and other particulate materials can exist in liquid and solid states [1]. The dominant paradigm for understanding liquid-solid transitions in these materials has been jamming [2][3][4][5][6], where the mechanical behavior is controlled by the ratio φ of total particle volume to total system volume. Jamming occurs when the system is isostatic, meaning that the number of contacts per particle Z is equal to the number required to constrain degrees of freedom in the system, Z = Z iso . At isostaticity, particles can no longer cooperatively rearrange to maintain a pressure p = 0, and further increasing φ requires compression of constituent particles such that p > 0 and Z > Z iso . Previous work [7][8][9][10][11] has shown evidence that jamming is a nonequilibrium critical transition, where a diverging cooperative length scale ξ J ∝ |φ − φ J | −νJ dominates the mechanical behavior near a critical volume fraction φ J , leading to Widom-like scaling relations [12] that relate p, φ − φ J , Z − Z iso , and other quantities.
However, jammed systems with p > 0 can also become liquid-like via yielding [13][14][15][16][17][18][19][20] if they are subjected to a shear stress τ such that Σ = τ /p is above the yield stress Σ c . Yielding in granular materials [21][22][23][24] and other jammed amorphous solids [25][26][27][28] also displays long-range cooperative effects: intermittent flow events can be triggered at Σ < Σ c if there is flow nearby, and Σ > Σ c can correspond to no flow if there is a nearby no-slip boundary [29]. Theories that include a diverging cooperative length scale ξ ∝ |Σ − Σ c | −ν can successfully predict both of these effects under time-averaged steady-state conditions. This suggests underlying critical behavior near Σ = Σ c , where local mechanically stable (MS) configurations of the material are correlated over ξ, but there is no particle-scale theory that explains the availability of MS states at varying Σ and over different length scales L. Moreover, since these studies are often carried out near φ = φ J , it is not known how yielding [30][31][32][33] relates to jamming or shear-induced jamming [34][35][36][37], which occurs near or below φ J .
In this Letter, we use numerical simulations to show that yielding is a distinct nonequilibrium phase transition from jamming. We study quasistatically sheared systems of repulsive, bidisperse disks and spheres at fixed dimensionless pressurep, wherep is varied over a wide range from near jamming to far above jamming. During shear, we measure the shear stress Σ, which increases during an initial shear regime and then plateaus, as stress is released in intermittent slips. As we vary the system size L, we find a scaling description at eachp with a diverging length scale ξ ∝ |Σ − Σ c (p)| −ν . Interestingly, we find two different values of ν: one during initial shear as stress builds up, ν ms ≈ 1.8 (in agreement with [38]), and another during the steady-state shear where stress is released in intermittent slips, ν slip ≈ 1.1 in 2D and ν slip ≈ 0.8 in 3D (consistent with, e.g., [33]). This scaling description, including the values of ν, is highly insensitive to the distance from jammingp. We find that Σ c (p) is also insensitive top forp < 10 −3 , but decreases logarithmically for higherp. 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 and spheres, 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 F ij = K(δ ij /|r ij | − 1)r ij , where δ ij is the average diameter of particles i and j, r ij is the vector displacement between the centers of particles i and j. Stresses are quantified by the Cauchy stress tensor, (1) where σ is the average Cauchy stress tensor, versus γ for a single simulation with 1,236 (orange) and 10,000 (black) particles. The inset shows a closeup of 0 ≤ γ ≤ 0.15. The first arrow indicates the initial shear strain γms required to find the first state at a particular value of Σ (in the case shown, Σ ≈ 0.077). Subsequent arrows denote the shear strains γ slip between states where the shear stress is less than the value of Σ being considered.
Here, α and λ are Cartesian coordinates, V is the system volume, r ij α is the α-component of the center-to-center separation vector between particles i and j, and F ij λ is the λ-component of the interparticle contact force. The sum over i and j includes all pairs of contacting particles.
Each simulation is prepared by placing disks randomly throughout the domain and then increasing the particle diameter D in small steps until reaching a target p = (σ xx + σ yy )/2. Using Lees-Edwards boundary conditions, we impose affine shear strain in small steps ∆γ = 10 −4 . At each shear step, the shear-periodic boundary is moved by ∆γ and y∆γ is added to the xposition 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 −Bv to each particle, where v is the absolute velocity of a given particle and B is the damping coefficient. We set B = 5 √ M p, 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 Σ = −2σ xy /(σ xx + σ yy ), as shown in Fig. 1(c). We measure Σ from 0 ≤ γ ≤ 3 in increments of ∆γ = 10 −4 for a total of 30,001 states per simulation. For each value of N andp, we simulate an ensemble of 400 systems.
We quantify distance above jamming byp = p/K, which gives an estimate of the relative overlap between particles (i.e.,p = 0.01 corresponds to particle-particle overlaps of roughly 0.01D). Figure 1(a) and (b)p = 0.2 andp = 0.001, respectively, with N = 24. Overcompression yields excess contacts such that ∆Z/Z iso > 0, where ∆Z ≡ (Z − Z iso ), Z is the number of contacts per particle, Z iso ≡ 2 (d (N − N r ) − d + 1) / (N − N r ) is the isostatic number of contacts, d is the number of spatial dimensions, and N r is the number of rattlers [4,[39][40][41]. Figure 2(a) shows ∆Z/Z iso plotted versus γ for N = 50 and varyingp. 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 the average value ∆Z /Z iso versusp for N = 50, 100, and 200, showing that the fraction of excess contacts is set bỹ p, nearly independent of system size [42].
Results.-As shown in Fig. 1(c), Σ increases with γ and then plateaus as potential energy is released in intermittent slips [30][31][32][43][44][45][46]. The fluctuations in Σ decrease with increasing N for allp, and we exploit the size scaling in these fluctuations (as in [33]) to demonstrate and quantify diverging spatial correlations near Σ c . Specifically, we use finite size scaling on three quantities related to the number density of mechanically stable configurations for a given Σ and N : (1) the cumulative distribution function F of states above a particular value of Σ during the slip avalanche regime, defined as γ > 0.5 (our results are insensitive to this choice); (2) the shear strain γ slip between MS states with an internal shear stress of at least Σ; and (3) the shear strain γ ms required to find the first MS state with an internal shear stress of at least Σ [38]. Figure 1(c) depicts γ ms and γ slip for a given Σ(γ) curve. Figure 3(a)-(c) shows these quantities plotted as functions of Σ or N , where ensemble averages are denoted with angle brackets. Figure 3 shows data forp = 0.05 in 2D, and is typical of allp in both 2D and 3D. As N is increased, the fluctuations in Σ decrease, and F approaches a step function, as shown in Fig. 3(a). Thus, MS states vanish sharply at some value Σ = Σ c (p) in the large system limit. Figure 3(b) shows γ slip , where we require at least one γ slip measurement per simulation. Our results are insensitive to this choice, unless the number of samples becomes very small. For Σ < Σ c , γ slip monotonically decreases with increasing N . For Σ > Σ c , γ slip first decreases and then increases with increasing N . Finally, γ slip is nearly independent of N for small Σ and increases strongly with N for larger Σ.
To collapse these curves, we posit a diverging length scale ξ ∝ |Σ − Σ c | −ν . In this case, finite size effects should enter through the quantity L/ξ, where L ≡ N 1/d with d being the number of spatial dimensions. An equivalent scaling can also be written using (Σ − Σ c ) L 1/ν ; see Refs. [9,38] for further discussion on similar systems. Figure 3(d)-(f) shows that the data in Fig. 3(a)-(c) col-lapses according to Here, f ms,± and f slip,± are dual-branch functions, with + and − denoting Σ above or below Σ c , respectively. Interestingly, we need two distinct values of ν for the initial strain, ν ms ≈ 1.8, and slip avalanche regime, ν slip ≈ 1.1 in 2D or ν slip ≈ 0.8 in 3D. The value ν slip ≈ 1.8 agrees with our previous result [38]. The difference between ν slip and ν ms suggests that there are important differences in how MS states are accessed between these two regimes. We obtain the critical parameters Σ c and ν slip by fitting the collapsed curves to appropriate functional forms using a Levenberg-Marquardt method [9,38]. We exclude system sizes with N < N min , and vary N min until our fits become insensitive to our choice of N min . We also performed the corrections-to-scaling analysis in [47], which yields the same result we find with the scaling forms in Eqs.
(2)-(4). For each value ofp, we first measure Σ c and ν slip by fitting f to the complimentary error function, f = 0.5 erfc (a(X − b)), where X = (Σ − Σ c ) L 1/ν slip , and Σ c , ν slip , a, and b are all fit parameters. A scaling collapse for F attained via this method is shown in Fig. 3(d), which uses 2D data at pressurep = 0.05. We apply similar methods to γ slip , shown in Fig. 3(b,e), by fitting data for γ slip to a low-order polynomial or by using crossingpoint analysis as in Ref. [48], both of which give similar results. We find similar values of Σ c and ν slip obtained from f , and we also find β slip /ν slip ≈ −1.05 for the case (2D,p = 0.05) shown in Fig. 3. Finally, we fit data for γ ms according to Eq. (4), and find with ν ms = 1.8 as in [38].
We then perform the same analysis for varyingp over the ranges 5·10 −5 <p < 2·10 −1 in 2D and 2·10 −4 <p < 2 · 10 −1 in 3D, spanning from near jamming (where ξ J is large) to far from jamming (where ξ J is small). The scaling description in Eqs. (2)-(4) and shown in Fig. 3 holds for all values ofp in both 2D and 3D. Figure 4 shows the critical parameters Σ c and ν slip plotted as a function of p. We find Σ c to be independent ofp forp < 10 −3 , and Σ c decreases logarithmically forp > 10 −3 . This occurs as excess contacts are added, which changes the structure of the force networks. Perhaps most importantly, the critical exponents during the slip avalanche regime are highly insensitive top, despite the wide variation in distance to jamming. Specifically, we find ν slip ≈ 1.1±0.1 in 2D, ν slip ≈ 0.8 ± 0.03 in 3D, β slip /ν slip ≈ −1 ± 0.1 in 2D, and β slip /ν slip ≈ −1.3 ± 0.1 in 3D. The uncertainty is estimated from the scatter in the data for differentp, as seen in Fig. 4. Other fit parameters, like a and b which enter into the fit for F as discussed above, are weakly dependent onp, especially forp < 10 −2 . We find a ≈ 1.2 and b ≈ −0.05 for low pressures, but forp > 10 −2 , a begins to increase while b does not change.
For the initial shear regime, we find that ν ms ≈ 1.8 is insensitive top. However, β ms /ν ms appears to vary from roughly 0.2 at highp to 0.6 at lowp. 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 [29,38,49] Discussion.-We have shown here that sheared amorphous soft sphere packings display finite size scaling that is consistent with a diverging length scale ξ ∝ |Σ−Σ c | −ν , and that this behavior is insensitive to distance from jammingp. Our results demonstrate that yielding is a distinct nonequilibrium critical transition from jamming. The value of Σ c varies asp is changed and extra contacts are added, but key features of the scaling functions, including the value of ν, are nearly independent of distance to jamming. Our results suggest yielding in, e.g., emulsions, foams, or granular materials may be less related to jamming as it is to yielding in other amorphous solids [30,32,33,43]. We note that ν slip ≈ 1.1 for 2D and ν slip ≈ 0.8 in 3D are similar to results from Ref. [33], which presented a scaling description for yield in amorphous materials and reported ν ≈ 1.1 for 2D and ν ≈ 0.7 for 3D. Figure 4(c) shows a phase diagram, similar to the τversus-φ jamming phase diagram from Refs. [2,34] and others, but withp on the horizontal axis and Σ = τ /p on the vertical axis. The solid blue line represents the critical yielding boundary in 2D from Fig. 4(a), and the solid vertical black line represents the jamming transition. Solid-like states exist only in the lower right region, denoted "jammed solid." Unjammed or fluid-like states can either hypostatic (Z < Z iso and p = 0) or hyperstatic (Z > Z iso and p > 0). We note that much previous work on critical scaling of jamming has studied yield stress behavior under shear at varied φ [8][9][10][11]. 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.  4(c). Future work may find other interesting behavior for very small pressures,p < 10 −5 , near the intersection of these two critical transitions.