{'fraction_non_alphanumeric': 0.044639408486453674, 'fraction_numerical': 0.02400264856811786, 'mean_word_length': 4.119774011299435, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dimensionality and viscosity exponent in shear-driven jamming

Collections of bidisperse frictionless particles at zero temperature in three dimensions are simulated with a shear-driven dynamics with the aim to compare with behavior in two dimensions. Contrary to the prevailing picture, and in contrast to results from isotropic jamming from compression or quench, we find that the critical exponents in three dimensions are different from those in two dimensions and conclude that shear-driven jamming in two and three dimensions belong to different universality classes.

Introduction A system of granular particles at zero temperature with contact-only interactions, undergoes a jamming transition, which is a transition from a liquid to a disordered solid, at a critical packing fraction φ J . As this is a phenomenon at zero temperature, there is no thermal equilibrium and it turns out that details of the jamming transition depend on the physical protocol by which the system jams; isotropic jamming and sheardriven jamming thus appear to be different phenomena.
Isotropic jamming results when the system is either compressed isotropically [1][2][3][4] or when it is rapidly quenched from T = ∞ to T = 0 at fixed volume [3,5,6]. In both cases the resulting jammed state has (in principle) an isotropic stress tensor. When compressed the particle packing φ is increased by slowly and isotropically compressing a system. As φ increases, particles come into contact with each other, at φ J a mechanically stable rigid backbone of particles percolates across the system, and the system jams. The precise value of φ J varies somewhat with the details of the protocol for compressing or quenching [3,7], as properties of the starting configurations and the rate of compression or quench.
In shear-driven jamming of frictionless particles the system is sheared at constant volume with a uniform shear strain rateγ. Below φ J -the jamming density of the shear-driven jamming transition, which is independent of the initial configuration-the system behaves as a liquid with a finite viscosity, limγ →0 (σ/γ), where σ is the shear stress. Above φ J a finite yield stress developes, limγ →0 σ > 0.
Early numerical simulations in 2D and 3D led to the conclusion that the critical exponents associated with isotropic jamming are independent of the dimensionality of the system [6]. More recently it has been demonstrated numerically that key non-trivial critical exponents for isotropic jamming agree quite well [8][9][10] with the values predicted analytically from an infinite-dimensional meanfield theory [11,12]. This observation has supported earlier claims that the upper critical dimension for isotropic jamming is d u = 2 [13,14], and that mean-field results apply for any d > d u . The prevailing view has been that the same should be true for shear-driven jamming [15], and theoretical models have been constructed that try to relate the critical exponents for shear-driven jamming to the mean-field values appropriate to isotropic jamming [15,16]. In this work we argue that this prevailing view is incorrect. By extensive numerical simulations, and a carefully quantitative analysis of the critical behavior, we show that the exponent associated with the diverging viscosity below φ J is clearly different in 2D and 3D, thus demonstrating that shear-driven jamming in physical dimensions cannot be considered a mean-field transition.
The expectation that shear-driven jamming in two and three dimensions should behave the same, seems to be taken over from the above-mentioned dimensionindependence found for isotropic jamming, together with the common result that weakly driven systems may be considered to be small perturbations about configurations in the absence of driving. This is however not applicable in the present situation since the shearing may never be considered to be a small perturbation. One way to see this is by considering the dimensionless friction at criticality which is µ ≡ σ/p ≈ 0.1 [17] (where p is pressure), which means that the system is highly anisotropic even in the limit of weak driving. A situation when linear response is applicable is in shearing simulations at finite temperatures and smallγ/T [18], but it is then found that linear response is applicable only as long as the system is close to isotropic, σ/p < 0.01.
Several attempts have been made to determine the critical behavior of shear-driven jamming [15,17,[19][20][21][22][23][24][25][26]. We here briefly review a few these methods. The first is to determine shear stress, σ(φ,γ) or pressure, p(φ,γ) from shear-driven simulations of soft disks at different densities and shear strain rates and make use of a scaling relation, described below [19,25], to try to extract the behavior in theγ → 0 limit. With the pressure-equivalent of the shear viscosity, η p ∼ p/γ, the divergence at the jamming density φ J is governed by the exponent β, (1) Since the particle overlaps get smaller for smallerγ, the limitγ → 0 is the hard particle limit. Another way to approach criticality is by doing shear-ing simulations with hard particles [17]. Since hard frictionless particles jam when the contact number is equal to z = z c = 2d (when z is determined after removing the rattlers from the system) the idea is to determine how η p diverges as z c is approached.
The advantage of this expression over Eq. (1) is that z c is known whereas φ J in Eq. (1) is unknown and has to be determined from the fitting together with the exponent. Eq. (2) therefore opens up a more direct analysis by just plotting η p vs δz ≡ z c − z.
For comparing determinations of β and β/u z one needs a value for u z , which in Ref. [21] was found to be u z = 1. This determination was however done without first removing the rattlers, and the precision has also been questioned [15]. Turning things the other way around, u z in 2D may be determined from β = 2.70 ± 0.15 from the scaling analysis [25,27] and β/u z = 2.69 ± 0.03 [28] (also shown in Fig. 1) which gives u z = 0.996 ± 0.057. Here and throughout the paper the quoted errors are max/min values, three standard deviations.
The essence of the shear-driven jamming transition is the slowing down of the dynamics, and the characterization of this dynamics is the idea behind a different but related method to study the jamming transition. In this method the ordinary shearing at a fixedγ is suddenly stopped and the system is made to relax to vanishing energy [28]. From the exponential relaxations one determines the relaxation time τ while one measures δz from the final configuration. (Note that this relaxation time is not the same as the relaxation time, commonly determined in steady state or at equilibrium, which is obtained from the self-part of the intermediate scattering function [29].) It turns out that τ determined from such relaxations behaves the same as η p [17,28] and we have The present paper presents shearing simulations of soft elastic particles and makes use of both the scaling analysis of pressure and the analysis of the relaxation time.
Models and simulations For the simulations we follow O'Hern et al. [6] and use a simple model of bi-disperse frictionless soft particles-disks or balls-in two and three dimensions with equal numbers of particles with two different radii in the ratio 1.4. Length is measured in units of the diameter of the small particles, d s . We use Lees-Edwards boundary conditions [30] to introduce a timedependent shear strain γ = tγ. We define the non-affine velocity, v i =ṙ i −v aff (r i ), obtained by subtracting off the uniform shear velocity v aff (r i ) ≡γy ix from the particle center of mass velocityṙ i . With r ij the distance between the centers of two particles and d ij the sum of their radii, the relative overlap is δ ij = 1 − r ij /d ij and the interaction between overlapping particles is V (r ij ) = ǫδ 2 ij /2; we take ǫ = 1. The force on particle i from particle j is The simulations are performed at zero temperature. We consider the interaction force f el i = j f el ij where the sum extends over all particles j in contact with i. The simulations discussed here have been done with the RD 0 (reservoir dissipation) model with the dissipating force We take k d = 1/2 and the time unit τ 0 = d 2 s k d /ǫ = 1/2. The equations of motion were integrated with the Heuns method with time step ∆t/τ 0 = 0.4. We simulate with N = 65536 particles and shear strain rates down tȯ γτ 0 = 10 −8 . By determining theγ below which finite size effects start become visible in additional simulations with N = 1024 and N = 4096, we conclude that our data with N = 65536 should not be affected by finite size effects.
Beside the ordinary simulations at constant shear strain rate we do two-step simulations where the constant shearing is suddenly stopped but the dynamics is continued, such that the systems relax to vanishing energy. From the exponential relaxations of p we determine the relaxation times τ and from the final configurations we determine the contact number z, after first removing the rattlers. The values of τ and δz ≡ z c − z from these relaxations will spread around averages that depend on both φ and the initialγ. It does however turn out when plotting points parametrically as τ vs δz, the points fall on a single common curve independent of the starting parameters [28]. This observation may be rationalized by considering that the final steps of the relaxation is probing the hard particle limit in which the dynamics is controlled by the properties of the contact network only, and thereby by the distance to jamming as measured by δz.
Results Our key result is summarized by Fig. 1 where τ vs δz in both two and three dimensions are shown by solid dots and open circles, respectively. The 2D data are from Ref. [28]. As always in the determination of critical exponents, we are interested in the data closest to criticality, i.e. at small δz; we note that the slopes at small δz in Fig. 1 are clearly different. Fitting data with δz < 0.08 to Eq. (2) gives the exponent β/u z = 3.35 ± 0.12 in 3D, clearly different from β/u z = 2.69±0.03 in 2D [28]. This is therefore strong evidence that shear-driven jamming in 3D and 2D belong to different universality classes. A more careful determination of the 3D exponent is given below.
Scaling analysis For a more detailed characterization of the critical behavior we turn to a scaling analysis of p(φ,γ). Following Ref. [25] the starting point is the scaling assumption below, where the second term gives the leading corrections-to-scaling, p(δφ,γ) = b −y/ν f (δφb 1/ν ,γb z ) + b −ω g(δφb 1/ν ,γb z ) . Here b is a length rescaling factor, f and g are scaling functions, ν is the correlation length exponent, z is the dynamical critical exponent, y is the scaling dimension of p, and ω is the correction-to-scaling exponent. Choosing b so thatγb z = 1, and with q = y/zν, this becomes p(δφ,γ) =γ q f p δφ γ 1/zν +γ ω/z g p δφ γ 1/zν .
We take f p and g p to be exponentials of sixth and third order polynomials in δφ/γ 1/zν , respectively. The data used for the fits are shown in Fig. 2 as η p (φ,γ) ≡ p/(γτ 0 ) for shear ratesγτ 0 = 10 −8 through 5 × 10 −5 . We start out by neglecting the corrections-to-scaling term and fitting to the simpler expression, We then adjust φ J , q, and zν together with the coefficients of the polynomial for f p , to get the best possible fit. As we don't know at the outset how big are the shear rates that can be used in the analysis, we do these fits with different ranges ofγ, takingγ min ≤γ ≤γ max withγ min τ 0 = 10 −8 andγ max τ 0 = 2.5 × 10 −7 through 5 × 10 −5 . The solid dots in Fig. 3 are from these fits. From the quality of the fits shown as χ 2 /DOF in panel (a) it is clear that the fits to the simpler Eq. (6) are good only when the data are restricted to very low shear rates. We then include corrections to scaling by fitting to the full expression, Eq. (5), taking ω/z and the coefficients of g p as additional free parameters. As this expression includes more fitting parameters we need more data in the fits and the analyses are therefore only done forγ max τ 0 ≥ 2.5 × 10 −6 . We conclude that the fit witḣ γ max τ 0 = 10 −5 gives reliable results by considering the quality of the fits together with the (weak) dependence oṅ γ max . We thus estimate β = (1−q)zν = 3.82±0.28, φ J = 0.6491 ± 0.0003, q = 0.233 ± 0.016, 1/zν = 0.200 ± 0.011, y = 1.16±0.03, and ω/z = 0.30±0.06. The quoted errors are max/min values (three standard deviations) whereas the error bars in the figures are ± one standard deviation. The errors are estimated with Jackknife resampling. The value β ≈ 3.8 in 3D is thus clearly different from the 2D value β ≈ 2.7 [25,27].
Corrections in the analysis of τ Due to the curvature of τ vs δz in 2D [28] it was found important to only make use of data for small δz in the determination of β/u z [28]. It was then found (not shown) that the determined β/u z increases as the range of δz decreases down to (δz) max = 0.08, but then stays stable. Decreasing (δz) max further only increases the statistical errors.
The analysis of the 3D data in Fig. 1 was similarly done by fitting to Eq. (3) with (δz) max = 0.08, and was indeed sufficient for demonstrating that this exponent is different in 3D compared to 2D. To check for the robustness of this determination, the lower left part of Fig. 4(a) shows β/u z vs (δz) max for the 3D data. In contrast to the behavior in 2D, this data does not clearly saturate but rather gives evidence for a trend to larger β/u z as (δz) max decreases. The value β/u z = 3.35 from Fig. 1 now only appears as a lower bound.
To try to get a better determination of the 3D exponent we now start from the assumption that the curvature in τ vs δz is related to corrections to scaling. By constructing a scaling expression for p/γ from Eq. (4), taking (−δφ) b 1/ν = 1 and δz ∼ (−δφ) uz , and noting that τ (φ,γ → 0) ∼ p(φ,γ → 0)/γ [28] one arrives at τ (δφ) = (δz) −β/uz f 0 + (δz) ων/uz g 0 , which is Eq. (3) with a correction term. Similarly to the scaling analysis which was done for differentγ max we fit our data with τ vs δz to Eq. (7) for δz ≤ (δz) max . As shown in Fig. 4(a), decreasing the range of data from (δz) max = 0.48 through 0.24 gives evidence for trends in both β/u z and ων/u z which appear to saturate at (δz) max = 0.30. We therefore read off β/u z = 3.7 ± 0.7. We also note that this value appears as a reasonable candidate to an extrapolation of the open squares in Fig. 4(a) to (δz) max = 0. Since we find numerically that β from p is equal to β/u z from τ , we conclude u z ≈ 1 in 3D, in agreement with the above-mentioned u z ≈ 1 in 2D.
Comparison with the literature Evidence for differing exponents in two and three dimensions has actually for some time been available in the literature. The first determinations of β/u z (there denoted by 1/δ) in Ref. [17], gave β/u z = 1/0.38 = 2.63 in 2D and β/u z = 1/0.34 = 2.94 in 3D. Those authors, however, did not consider this a significant difference. The main source of uncertainty in these analyses is whether the data are sufficiently close to criticality to give the true critical behavior. In a later paper by the same group [15], simulations closer to criticality-i.e. at smaller δz-gave β/u z = 1/0.3 ≈ 3.3 in 3D, but as that paper was focused on comparisons with theory they didn't comment on possible differences between two and three dimensions. Their values do however agree nicely with our analyses in Fig. 1.
We also note from Fig. 3(b) that a simple scaling analysis of p(φ,γ), based on Eq. (6), without corrections to scaling, gives β ≈ 2.9, close to the 2D value, provided that one includes in the fit larger values of the strain ratė γ. The bad quality of the fit in Fig. 3(a) makes clear that these low values cannot be correct but it shows that analyses may give the erroneous conclusion that the critical behavior in 3D and 2D are the same, seemingly confirming the prevailing paradigm. Another example of a low value in the literature is β ≈ 1/0.391 = 2.56 [26]. This, again, appears to be an effect of using data too far from criticality, as the fits, according to their Fig. 4(a), include points for densities down to, or below, φ ≈ φ J − 0.05. To compare, the data used in the scaling analyses in the present work are restricted to |φ − φ J | ≤ 0.017.
Discussion Recent attempts by the group of Wyart to determine the exponents analytically, in terms of the exponent θ e , rely on examining the properties of the opening and closing of contacts [15,16]. The exponent θ e characterizes the distribution of weak forces in packings from isotropic jamming and has been found to be θ e ≈ 0.42311 by analytic calculations in infinite dimensions [11,12]. It is also found to be the same in 2D and 3D [8][9][10] and is believed [15] to be the same also in the shear-driven case.
The result of the present Letter, that critical exponents for shear-driven jamming are different for 3D compared to 2D, is however in conflict with a picture where the exponents only depend on the dimension-independent exponent θ e . One possible reason for this difference could be that their relations [15,16] describe the typical particle motion whereas the dissipation (and the viscosity) is instead dominated by a small fraction of particles with the highest velocity-a fraction which decreases as jamming is approached [32]. The investigation into this issue appears as an important direction for future research.
Conclusion From shear-driven simulations of elastic particles in three dimensions together with previous results for two dimensions, we determine the critical exponents of the shear-driven jamming of frictionless athermal particles and conclude-in variance with the prevailing picture-that the 3D and 2D transitions do not belong to the same universality class. I thank S. Teitel for suggestions, discussions, and a critical reading of the manuscript. Simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at HPC2N.