J ul 2 01 9 Critical scaling for yield is independent from distance to isostaticity

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 μ of these states obey a scaling description with diverging length scale ξ ∝ |μ−μc| −ν . 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.

Prior studies on soft sphere packings, which are commonly used to model these materials, have framed longrange cooperative behavior in terms of a nonequilibrium critical transition that occurs at the isostatic point, also called "point J" [2,[22][23][24][25][26][27][28].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 iso .This occurs at a given volume fraction φ = φ c in the large-system limit.At isostaticity, p = 0, but further compression leads to increasing p.A cooperative length scale ξ J ∝ |φ − φ J | −νJ diverges at the isostatic point, which then controls the mechanical response and leads to Widom-like scaling relations [29] that relate p, (φ − φ J ), (Z − Z iso ), and other quantities.ξ J is large near isostaticity (i.e., small p), characterized by an excess of spatially extended, low-energy modes of the system [30,31].For increasing p, ξ J decreases, leading to more localized modes as well as smaller and more localized particle rearrangements.
In contrast, nonlocal rheological descriptions of jammed materials [16,[18][19][20][21]32] often include a diverging cooperative length scale that depends not on packing fraction but on distance to a critical shear stress, i.e., ξ ∝ |µ − µ c | −ν .These rheological models, including our previous paper [32], describe materials that are near φ = φ J , so it is not known how the cooperative length scale underlying these models relates to the iso-static 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 p fixed and measuring µ, which increases during an initial shear regime and then plateaus as stress is released in intermittent slips.The statistics of µ obey a scaling description with a diverging length scale ξ ∝ |µ−µ c | −ν , where ν ms ≈ 1.8 during initial shear buildup (in agreement with [32]) and ν slip ≈ 1.1 in two dimensions (2D) and ν slip ≈ 0.8 in three dimensions (3D) during the intermittent slip regime.The scaling functions and the values of ν are highly insensitive to the distance from isostaticity set by p, which we vary over nearly four orders of magnitude.µ c (p) is insensitive to p for p ≤ 10 −3 , but decreases logarithmically for higher 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 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, 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 particles 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 µ = −σ xy /p, 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 and p, we simulate an ensemble of 400 systems.We quantify distance above isostaticity by p = 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 Figure 2(a) shows ∆Z/Z iso plotted versus γ for N = 50 and varying 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 ∆Z /Z iso versus p is similar for N = 50, 100, and 200.Thus, the fraction of excess contacts and thus the distance to isostaticity is set by p, nearly independent of system size [36] or the presence of shear deformation [37].Scaling near yield.-Asshown in Fig. 1(c), µ increases with γ and then plateaus as potential energy is released in intermittent slips [38][39][40][41][42][43][44].This curve represents a series of jammed states that the system passes through while sheared.The fluctuations in µ decrease with increasing N for all p, and we exploit the size scaling in these fluctuations (as in [45]) to demonstrate and quantify diverging spatial correlations near µ 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 µ and N : (1) the cumulative distribution function F of states above a particular value of µ during the slip avalanche regime, defined as γ > 0.  (our results are insensitive to this choice); (2) the shear strain γ slip between mechanically stable (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 µ. 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.The data shown in Fig. 3 represents only a single value of p = 0.05 in 2D, but it is typical of all p in both 2D and 3D, as we demonstrate below in Figs. 4 and 5.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.[27,32] for further discussion on similar systems.Figure 3(d)-(f) shows that the data in Fig. 3(a)-(c) collapses 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 [32], which was only calculated near to isostaticity; in Fig. 3, it is calculated far from isostaticity.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 [27,32].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 [46], which yields the same result we find with the scaling forms in Eqs. ( 2)-(4).
We then perform the same analysis for varying p 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 isostaticity (where ξ J is large) to far from isostaticity (where ξ J is small).The scaling description in Eqs. ( 2)-( 4) and shown in Fig. 3 holds for all values of p in both 2D and 3D.We show data for an additional pressure in 2D, p = 0.0005, in Figure 4. We also show data in 3D with 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 µ c and ν slip plotted as a function of 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 ≤ N ≤ 1024) with 400 simulations per system size, so the plateau in Fig. 6 is not a system size effect.We find µ c to be independent of p for p ≤ 10 −3 , and µ c decreases logarithmically for p > 10 −3 , which agrees with Favier de Coulomb et al. [37].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 The uncertainty is estimated from the scatter in the data for different p, as seen in Fig. 6.For the initial shear regime, we find that ν ms ≈ 1.8 is insensitive to p.However, β ms /ν ms appears to vary from roughly 0.2 at high p to 0.6 at low 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 [32,47,48].
Discussion.-We have shown here that sheared amorphous soft sphere packings display finite size scaling that is consistent with a diverging length scale ξ ∝ |µ − µ c | −ν .The value of µ c varies as 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 p. Considering the correlation length for isostaticity ξ J ∝ |φ − φ c | −νJ , if one assumes that ν J is order unity [27] and p ∝ (φ − φ c ) for harmonic interactions [2], then varying p over this range represents ξ 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 ν slip ≈ 1.1 for 2D and ν slip ≈ 0.8 in 3D are similar to the values ν ≈ 1.1 for 2D and ν ≈ 0.7 for 3D from Ref. [45], which presented a scaling description for yielding in amorphous materials [38,40,43].tical 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 iso and p = 0) or hyperstatic (Z > Z iso and p > 0).Some previous work on critical scaling near isostaticity has studied the onset of yield stress behavior under shear at varied φ [24,25,27,28,49].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.

FIG. 1 .
FIG. 1. (a, b) Illustrative snapshots during Lees-Edwards shear at γ = 0.33 with dimensionless pressure (a) p ≡ p/K = 0.2 (far from isostaticity) and (b) p = 0.001 (near isostaticity).(c) Plot of dimensionless shear stress µ 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.

FIG. 2 .
FIG. 2. (a) Excess contacts ∆Z/Ziso versus the shear strain γ for individual simulations with N = 50, where ∆Z ≡ Z − Ziso, Z is the coordination number once rattlers have been removed, and Ziso is the coordination number for an isostatic system.(b) Mean ∆Z/Ziso versus p over 50 simulations for N = 50, 100 and 200, showing that p gives the fraction of extra contacts, independent of system size.
(a) and (b) show p = 0.2 and p = 0.001, respectively, with N = 24.Overcompression yields excess contacts such that ∆Z/Z iso ∼ √ p, 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, 33-35].

FIG. 4 .FIG. 5 .
FIG. 4. Scaling collapses in 2D at low dimensionless pressure p = 0.0005, 32 ≤ N ≤ 512.Fraction of states F above µ, unscaled (a) and scaled (b), with µc = 0.097 and ν slip = 1.10.γ slip versus N , unscaled (c) and scaled (d), with µc = 0.097, ν slip = 1.10, and β slip /ν slip = −1.05.In both cases, the unscaled plots show selected values of N or µ while the scaled plots show all data.(a) (b) Figure 4. We also show data in 3D with p = 0.05 in Fig.5.In both cases, the scaling functions are almost indistinguishable from those shown in Fig.3.Figure6shows the critical parameters µ c and ν slip plotted as a function of p.Each data point in Fig.6represents a fit of all data (as in Figs. 3 through 5) over many system sizes (typically 32 ≤ N ≤ 1024) with 400 simulations per system size, so the plateau in Fig.6is not a system size effect.We find µ c to be independent of p for p ≤ 10 −3 , and µ c decreases logarithmically for p > 10 −3 , which agrees with Favier de Coulomb et al.[37].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 p, despite the wide variation in distance to isostaticity.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 different p, as seen in Fig.6.For the initial shear regime, we find that ν ms ≈ 1.8 is insensitive to p.However, β ms /ν ms appears to vary from roughly 0.2 at high p to 0.6 at low 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[32,47,48].Discussion.-Wehave shown here that sheared amorphous soft sphere packings display finite size scaling that is consistent with a diverging length scale ξ ∝ |µ − µ c | −ν .The value of µ c varies as p is changed and extra contacts are added, but the forms of the scaling functions (as shown in Figures3 , 4, and 5) and the values of the critical exponents are nearly independent of distance to isostaticity over nearly four orders of magnitude in p. Considering the correlation length for isostaticity ξ J ∝ |φ − φ c | −νJ , if one assumes that ν J is order unity[27] and p ∝ (φ − φ c ) for harmonic interactions[2], then varying p over this range represents ξ 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 ν slip ≈ 1.1 for 2D and ν slip ≈ 0.8 in 3D are similar to the values ν ≈ 1.1 for 2D and ν ≈ 0.7 for 3D from Ref.[45], which presented a scaling description for yielding in amorphous materials[38,40,43].Figure6(c) shows the Liu-Nagel jamming phase diagram from, e.g., Refs.[2,49,50] and many others, but with p on the horizontal axis and µ = τ /p on the ver-

Figure 6 (FIG. 6 .
Figure 4. We also show data in 3D with p = 0.05 in Fig.5.In both cases, the scaling functions are almost indistinguishable from those shown in Fig.3.Figure6shows the critical parameters µ c and ν slip plotted as a function of p.Each data point in Fig.6represents a fit of all data (as in Figs. 3 through 5) over many system sizes (typically 32 ≤ N ≤ 1024) with 400 simulations per system size, so the plateau in Fig.6is not a system size effect.We find µ c to be independent of p for p ≤ 10 −3 , and µ c decreases logarithmically for p > 10 −3 , which agrees with Favier de Coulomb et al.[37].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 p, despite the wide variation in distance to isostaticity.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 different p, as seen in Fig.6.For the initial shear regime, we find that ν ms ≈ 1.8 is insensitive to p.However, β ms /ν ms appears to vary from roughly 0.2 at high p to 0.6 at low 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[32,47,48].Discussion.-Wehave shown here that sheared amorphous soft sphere packings display finite size scaling that is consistent with a diverging length scale ξ ∝ |µ − µ c | −ν .The value of µ c varies as p is changed and extra contacts are added, but the forms of the scaling functions (as shown in Figures3 , 4, and 5) and the values of the critical exponents are nearly independent of distance to isostaticity over nearly four orders of magnitude in p. Considering the correlation length for isostaticity ξ J ∝ |φ − φ c | −νJ , if one assumes that ν J is order unity[27] and p ∝ (φ − φ c ) for harmonic interactions[2], then varying p over this range represents ξ 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 ν slip ≈ 1.1 for 2D and ν slip ≈ 0.8 in 3D are similar to the values ν ≈ 1.1 for 2D and ν ≈ 0.7 for 3D from Ref.[45], which presented a scaling description for yielding in amorphous materials[38,40,43].Figure6(c) shows the Liu-Nagel jamming phase diagram from, e.g., Refs.[2,49,50] and many others, but with p on the horizontal axis and µ = τ /p on the ver-