Characterizing the Nature of the Yielding Transition
aa r X i v : . [ c ond - m a t . s o f t ] O c t Characterizing the Nature of the Yielding Transition
S. H. E. Rahbari , ∗ J. Vollmer , and Hyunggyu Park School of Physics, Korea Institute for Advanced Study, Seoul 130-722, South Korea and Faculty of Mathematics, Georg-August Univ. Gttingen, 37077 Gttingen, Germany
Particulate matter, such as foams, emulsions, and granular materials, attain rigidity in a denseregime: the rigid phase can yield when a threshold force is applied. The rigidity transition inparticulate matter exhibits bona fide scaling behavior near the transition point. However, aprecise determination of exponents describing the rigidity transition has raised much controversy.Here, we pinpoint the causes of the controversies. We then establish a conceptual framework toquantify the critical nature of the yielding transition. Our results demonstrate that there is aspectrum of possible values for the critical exponents for which, without a robust framework, onecannot distinguish the genuine values of the exponents. Our approach is two-fold: ( i ) a precisedetermination of the transition density using rheological measurements and ( ii ) a matching rulethat selects the critical exponents and rules out all other possibilities from the spectrum. Thisenables us to determine exponents with unprecedented accuracy and resolve the long-standingcontroversy over exponents of jamming. The generality of the approach paves the way to quantifythe critical nature of many other types of rheological phase transitions such as those in oscillatoryshearing. I. INTRODUCTION
Yield stress materials such as toothpaste, hair gel,mayonnaise, and cement, are ubiquitous. These mate-rials are used in pharmaceutical and cosmetics manu-facturing, as well as the oil, concrete, and food indus-tries [1]. Because of their wide applicability in everydaylife, a quantitative description of their rheological behav-ior is pivotal. The physical origin of the yield stress de-pends on the microscopic details of the system and can beclassified into three main categories: dynamic arrest inBrownian suspensions known as the glass transition [2],mechanical (meta)stability in athermal systems or jam-ming [3, 4], and attractive interactions [5, 6]. Thixotropicyield stress fluids [1], which exhibit memory effects anda bifurcation in the viscosity, are outside the focus of thecurrent study.The relation between shear stress σ and shear rate ˙ γ ,known as a flow curve in a yield stress material, can bedescribed as a Herschel-Bulkley (HB) relation: σ = σ y + K ˙ γ ∆ , (1)where σ y is the yield-stress and ∆ < η = σ/ ˙ γ . As a result of the threshold σ y ,the viscosity of a yield stress material diverges for ˙ γ → et al. [8]repeated the same experiment and showed that thosemeasurements at low shear stresses never reached a ∗ [email protected] stationary state and that the apparent finite viscositywas an artifact of the measurement.A consensus regarding the existence of the yieldstress has emerged. However, the technical difficultiesof its measurements remain a challenge. Despite theseadvances in the understanding of the yield stress, adescription of the non-linear flow curves in the fluidstate remains an open problem. In the traditionalapproach, the shear-thinning exponent ∆ is obtainedby a power law fit to σ − σ y versus ˙ γ . However,recent numerical simulations showed that σ − σ y vs ˙ γ exhibits two distinct scaling regimes described by twodifferent exponents, ∆ and ∆ ′ , for small and large shearrates [9–11], respectively. As a result, fitting a HB-typerelation to such data will be prone to pitfalls due toa bias towards larger shear rates, which in turn willgive rise to a misleading quantification of the flow curves.The problem becomes even more dramatic for the caseof matter with granularity [12]. Soft particulate mat-ter, such as gels and emulsions, flow freely in the diluteregime and attain yield stress above a threshold density φ J in the dense regime. This yielding transition exhibitsa rich class of scaling behavior of the flow curves de-scribed by critical exponents (we will give a brief overviewabout different scaling regimes of the rigidity transitionin the next section). Despite many efforts by differentgroups [13–18], a precise determination of the critical ex-ponents remains disputed.Here, we establish a conceptual framework for the scal-ing quantification of the flow curves of a wide range ofyield stress materials. We resolve the long-standing dis-pute over exponents of the rigidity transition. II. RIGIDITY TRANSITION: A BIRD’S EYEVIEW
Depending on the shear rate and packing fraction,soft frictionless spheres display a rich phenomenology ofdistinct rheological regimes. This makes soft frictionlessspheres
Drosophila of particulate matter.In the dilute regime of particulate materials, flowcurves at small shear rates are given by σ ∝ ˙ γ n where n = 1 and 2 for Newtonian and Bagnoldian scalings, re-spectively. Because soft particles barely deform at smallshear rates, σ ∝ ˙ γ n corresponds to the so-called hard-corelimit. The exponent n has been shown to depend on theReynolds number of the system such that for overdampedsystems the Newtonian regime (non-inertial) is recoveredand for n = 2 the system must be under-damped (iner-tial) [19]. The transport coefficient, which is given byshear viscosity η = σ/ ˙ γ n , at ˙ γ →
0, depends only onthe packing fraction φ and diverges upon approachingthe jamming density η ∝ | δφ | − β , where δφ = φ − φ J isthe distance from jamming. The exponent β also char-acterizes the hard-core limit of the system. Accordingly,this exponent must be independent from the microscopicdetails of the system [20, 21]. At φ = φ J , the system ex-hibits pure power-law rheology σ ∝ ˙ γ q with q < φ > φ J , the system displays threshold rheology and flowcurves that may be described by the HB model givenby Eq. 1. In this model, the shear-thinning exponent ∆is shown to be related to the behavior of the system inthe hard-core limit at φ < φ J and thus to the exponent β [20]. The yield stress also scales with the distance fromjamming σ y ∝ δφ y .As one can see, upon approaching the jamming point,the rheology changes dramatically due to the collectivebehavior of particles [22]. Consequently, the rheology canno longer be described by trivial exponents such as n = 1or 2 and thus the system becomes shear-thinning with anon-trivial scaling dimension q <
1. This is a signatureof a growing length scale in the system [9, 23, 24], whichis the hallmark of critical phenomena. Even thoughthis system is non-equilibrium and athermal, Olsson andTeitel [25] used renormalization group formalism [26] ofequilibrium phase transitions to capture the critical na-ture of this dynamic transition. The jamming point at δφ = 0 , ˙ γ → , T = 0 , and L → ∞ is a genuine dynamiccritical point.Altogether, any of the above scaling limits can be re-trieved by choosing appropriate limits of a scaling func-tion F and an arbitrary length scale b in the followingscaling ansatz: (derivation given in Appendix B): σ ( δφ, ˙ γ, L, w ) = b − y/ν F ( δφ b /ν , ˙ γ b z , L − b, w b − ω ) , (2)where F is a homogeneous scaling function, L is thesystem size, and w is an auxiliary variable. This scal-ing ansatz is traditionally used to find relations between authors y q ω/z Otsuki & Hayakawa (theory) [14] 1 2 / − Hatano [13] 1 . . − Hatano [15] 1 . . − Otsuki & Hayakawa (simul.) [16] 1 .
09 0 . − DeGiuli et al. [17] 1 0 . . et al. [18] 1 . . . et al. [28] 1 − − TABLE I. Critical exponents reported by different authors.As the data get closer to the critical point, exponents system-atically change (for a comprehensive discussion, cf. [18]). different exponents. Inserting b = ˙ γ − /z for L → ∞ inEq. 2, we arrive at: σ = ˙ γ q F (cid:18) δφ ˙ γ q/y (cid:19) , (3)where q = y/ ( zν ). Here, we assume proximity of thecritical point where the auxiliary variable w can be ne-glected.The immediate outcome of Eq. 3 is that all the datamust collapse into a master curve when plotted σ/ ˙ γ q vs δφ/ ˙ γ q/y , providing that three free parameters q , y ,and φ J , are fine tuned. Notably, in the early stage ofthis topic, this method, ı.e., collapse of the data, hasbeen extensively used by many authors to estimate q , y ,and φ J [13, 15, 16, 25, 27]. A summary of the existingpredictions for these exponents is given in Tab. I.These reports were not conclusive because of the largerange of reported exponents and critical densities. Thereason for this was because the quality of the collapseswere judged based on the visual appeal of the plots.Later, Olsson and Teitel used a quantitative method tocompute the quality of the collapses. The method wasbased on ( i ) exponential parametrization of the scalingfunction F ( x ) = exp (cid:16)P n =0 a n x n (cid:17) and ( ii ) going intounprecedented small shear rates down to 10 − in thedimensionless scale [9, 18]. However, the expansion of F ( x ) may be prevented because, as x → F ( x ) maynot be analytic. Also, for reasons that we describe inthe next paragraph, going into shear rates as small as˙ γ ≈ − contaminates the scaling behavior.It is well known that in the jammed state a shearedparticulate system exhibits shear localizations, alsoknown as shear-transformation-zones [29, 30]. Thesestress anomalies relax through long-range system-wideavalanches. Each avalanche can trigger other activezones that will in turn result in a domino of plasticevents and relaxations. At very small shear rates, theseavalanches are globally correlated and poise the systeminto an effective critical state [31, 32]. This resultsin scale-free distributions of avalanches with exponentsthat are generally smaller than 2 [33, 34]. To obtaina flow curve σ ( ˙ γ ), one should perform time averagingfor shear stress over the time series. However, dueto the scale-free distribution of avalanches with theaforementioned range of exponents, the first and secondmoments of the shear stress cannot be well-defined.Consequently, the time averaged shear stress at verysmall shear rates possesses error bars that are as largeas the average values.To avoid the above problems, we describe a generalframework that requires neither data collapse nor expan-sion of the scaling function. Additionally, measurementsare performed outside the avalanche region, ı.e., notat very small shear rates. These simplify the problemdramatically and enable us to resolve the controversyover exponents. Our approach is two-fold: first, inSec. III, we describe how we nail down the criticaldensity. Second, in Sec. IV, we present our matchingrule that selects critical exponents from a wide spectrumof possible values. III. HUNT FOR φ J Precise determination of the critical exponentsstrongly depends on whether the critical density φ J isaccurately determined. In this section, we explain howwe nail down the transition density φ J using rheologicaldata. To achieve this goal, we define successive slopes ofthe flow curves m as: m = d ln σd ln ˙ γ , (4)where d stands for the derivative. This can be easilycalculated from Eq. 3: m = q − qy δφ ˙ γ q/y F ′ ( x ) F ( x ) , (5)where x = δφ/ ˙ γ q/y , F ′ ( x ) = d F ( x ) /dx .Eq. 5 provides an immediate prediction: if one plots m vs ˙ γ for different packing fractions, exactly at jammingdensity δφ = 0, the successive slope for all shear rateswill be equal to the critical shear-thinning exponent m = q . For δφ >
0, the successive slope convergesto m = q at large shear rates and deviates fromthat value for ˙ γ → γ − q/y . Similar be-havior is expected for δφ < φ J : thecritical density is given by a horizontal line of the m − ˙ γ dependence that distinguishes off-critical densitieswith opposite curvatures. However, it is practicallyimpossible to recover a straight horizontal line for m at φ J in the critical region of ˙ γ →
0. This is due toelasto-plastic critical fluctuations near the critical point,which we mentioned in Sec. II.The remedy for this problem is to stay away from theregion where the successive slope displays huge fluctu-ations. In such a regime, correction-to-scaling must betaken into account. From Eq. 2, the leading correction-to-scaling term at φ = φ J reads: σ = ˙ γ q (cid:16) c + c ˙ γ ω/z (cid:17) , (6)where c and c are constants and ω/z is the lead-ing correction-to-scaling exponent (see Appendix B forderivation). For off-critical densities φ = φ J , an extraterm proportional to δφ must be added to Eq. 6. Thisterm again has an inverse algebraic dependence on ˙ γ sim-ilar to that in Eq. 5. One can easily calculate the corre-sponding successive slope of Eq. 6 as: m = q + k ˙ γ ω/z , (7)where q is the asymptotic exponent and k is a constant.This shows the behavior of the successive slopes at φ J ,which distinguishes that of off-critical densities withopposite curvatures.Now let us calculate the asymptotic values of the suc-cessive slopes for different densities at ˙ γ →
0. For δφ < σ ∝ ˙ γ n , which results in m = n . At δφ = 0, σ ∝ ˙ γ q , then m = q . For δφ >
0, the yield stress emerges, whichamounts to a dependence ∝ ˙ γ and thus m = 0. Insummary: lim ˙ γ → m = n φ < φ J q φ = φ J φ > φ J . (8)We summarize the behavior of the successive slope offlow curves in a schematic diagram in Fig. 1. This dia-gram demonstrates the simplicity behind our frameworkto find φ J . In a semi-log plot of m vs ˙ γ , all of the sub-and super-critical densities curve in opposite directions,except at φ J .In our strategy to find φ J , we first obtain flow curvesfor an intermediate system size. We mark the range ofdensities where the curvature of the successive slopeschanges. We then zoom into the region by simulating alarger system size and nail down φ J . Finally, we checkwhether our estimated φ J is robust against finite-sizeeffects.We perform extensive large-scale two dimensionalmolecular dynamics simulations of frictionless disks in asimple shear flow. In our simulations, we dissipate thenormal component of the relative velocity of collidingparticles. This dissipation law leads to Bagnoldianscaling in the dilute regime. The Newtonian regime isrecovered when the transverse component of the relativevelocity is dissipated [21]. This regime is not explored inthis work. Further details of the simulations are givenin Appendix A. In Fig. 2, we display the successiveslope m vs shear rate ˙ γ for different packing fractions φ for a system of intermediate size L = 100. Thecurvature of the curves changes in the range between φ = 0 .
843 and 0 . φ J . We will zoom into this region to determine φ J with a higher resolution and larger system sizes. All ofthe curves corresponding to different packing fractionsshow a tendency to converge at large shear rates. Thisis in accord with the predictions by Eq. 5. One cansee that upon decreasing the shear rate, the far-topcurves show a tendency to converge towards the valueof the asymptotic exponent n = 2 and the far-bottomcurves to 0. This is again in accord with the predictionby Eq. 8. A dashed line shows an estimation for thevalue of q = 0 .
6. We note that this line tends towardssmaller values upon increasing L . For L ≥ q does not change. We note thatfor ˙ γ < − , the successive slope in the critical rangeof densities displays giant fluctuations reminiscent ofcritical fluctuations. We observe these fluctuations forsystems of larger spatial extents for ˙ γ < − . Therefore,in the rest of the paper, we do not consider data with˙ γ < − in our analysis. To summarize the results for L = 100, the crude estimation for the transition densityis φ J ≈ . ± . m q0 0 FIG. 1.
Schematics of successive slope.
A schematicfigure shows the successive slope vs shear rate in a semi-logscale. According to Eq. 5, curves at sub- and super-criticaldensities have opposite curvatures. In the asymptotic limit˙ γ →
0, all curves converge to the asymptotic exponents givenby Eq. 8. The curve corresponding to φ J is the only curve that does not bend upward or downward and whose offset isequal to the non-trivial critical exponent q . critical exponent is q ≈ .
6. Next, we will zoom into thecritical region with substantially larger system sizes tofind φ J .According to elasticity theory, shear stress and pres-sure are both components of a single entity known as thestress tensor. Different components of the stress tensorprovide information about momentum transfer in differ-ent directions into/along imaginary surfaces in the sys-tem [35]. However, whether the shear stress and pressurescale equivalently with shear rate is not at all an obviousfact. According to Peynneau and Roux [36] and morerecently by Baity et al. [37] a finite stress anisotropy, δp ∝ p xx − p yy , gives rise to a small rotation of princi-pal axes of of the stress tensor from those given by thestrain tensor. This gives rise to distinct scaling of theshear stress and pressure when there is a stress anisotropyin the system; this usually happens at high shear rates.However, the stress anisotropy is negligible at small shearrates near jamming [38]. This is also confirmed by ourresults. Thus, it is a widely accepted fact that the asymp-totic scaling of the shear stress and pressure are equiva-lent. This assumption has been adopted by many recentstudies, cf. [18]. More recently, Suzuki and Hayakawaprovided a rigorous derivation of this based on a µ - J -6 -5 -4 m = d l n σ xy / d l n γ ˙ shear rate FIG. 2.
Successive slope.
Successive slope m vs shear rate˙ γ for different φ . The value of each φ is given in the legend.The number of particles varies from N = 7183 to 7742 for φ = 0 .
835 and 0 .
90, respectively. The spatial extent of thesystem is L = 100. The curvature of the curves changes inthe range of φ = 0 .
843 and 0 . φ J and the range for refined measurements witha much larger resolution on φ for substantially larger systemsizes. rheology [39]. We will use this assumption in the nextsection to nail down the critical exponent q .We display refined measurements in Fig. 3 for differentsystem sizes up to L = 300. Panel a and b refer tothe successive slopes of the shear stress and pressure,respectively. One can see that for all densities, there is astrong system size dependence for L < L ≥ φ = 0 .
843 and 0 .
844 clearlyhave opposite curvatures for all system sizes. We zoominto this region to find the critical density. Filled squarescorrespond to φ = 0 . L = 300. These dataare averaged over 7 different ensembles. The rest of thedata are obtained from a single realization. For L = 300,a closer inspection of data at φ = 0 . . φ = 0 . φ = 0 . φ = 0 . φ J = 0 . ± . et al. [40] and Vagberg et al. [18].A closer inspection of the successive slope of shearstress σ (panel a ) and pressure p (panel b ) reveals astronger corrections-to-scaling of the shear stress. Here,stronger corrections-to-scaling means a larger amplitudeof the scaling function of Eq. 7. However, as we havementioned in the previous paragraph the asymptotic ex-ponents must be equivalent for both pressure and shearstress. Interestingly, a stronger corrections-to-scalingof shear stress has been reported by other authors [9, 18].One can see that φ J does not have a strong dependenceon the system size. However, the asymptotic exponentchanges continuously from 0 . . L = 50 to 300, respec-tively. Estimation of q for L = 300 is not straightforwardbecause of the complexity of the scaling function for largesystem sizes, ı.e., the dependence of m to ˙ γ . In the nextsection, we describe a systematic method to nail downthe critical exponents. IV. HUNT FOR EXPONENTS
In this section, we describe how we nail down thecritical exponents. The easiest way to find criticalexponents is to obtain them via fitting Eq. 7 to thesuccessive slope curve at φ J in Fig. 3 using q , k , and ω/z as free fitting parameters. We call this a blind fitting. Notably, a 3-parameter fitting correspondsto the optimization of a residual function in a 3 + 1dimensional space. This function is rugged and hasmany local basins. Each fitting algorithm/software willfind one such local minimum. This will cause a zoo ofdifferent values for the exponents due to the rugged nature of the residual function.To avoid fitting artifacts, we hold the correction-to-scaling exponent ω/z fixed and obtain the asymptoticexponent q via fitting a linear function to m versus ˙ γ ω/z .We vary ω/z in a range between 0 . .
5, and werecord the corresponding q . The contour lines of the fitsare given in Fig 4-a for both shear stress and pressure.Each contour line represents all possible outcomes of q and ω/z via a 3-parameter blind fitting. Each pointon the contour lines corresponds to one basin. Now,the crucial question becomes about which point on thecontour line can be considered the corresponding pointfor correct exponents.As we previously noted, pressure and shear stresscan have different scaling functions; however, asymptoticcritical exponents must strictly be equivalent. This pro-vides a matching rule , which allows us to pick up thecorrect exponents based on the crossing point of the con-tour lines of p and σ . Fig. 4-a demonstrates that such amatching point really does exist, and we read exponents q = 0 .
41 and ω/z = 0 .
365 for both p and σ . We notethat within error bars, the crossing point gives the same q for L ≥ ω/z is not stable. Therefore, weperform the finite-size scaling analysis for ω/z with fixed q ( L = ∞ ) = 0 .
41 via m ( L ) − q ( L = ∞ ) ∝ ˙ γ ω/z ( L ) . (9)We fit Eq. 9 and obtain ω/z as a function of L . Weplot ω/z ( L ) vs L − in Fig. 4-b. One can see that ω/z levels off at 0 .
35 for the largest system sizes. This givesus the asymptotic value of the leading correction-to-scaling exponent ω/z ( L = ∞ ) = 0 .
35. We summarizedthe values of critical exponents in Tab. II.Having obtained both q ( L = ∞ ) and ω/z ( L = ∞ ), wearrive at our final vital benchmark. We now hold q and ω/z fixed to their asymptotic L = ∞ values and fit Eq. 7to the data to obtain the amplitude k . The resultingcurves are shown as solid lines in both Fig. 3-a and -b.We obtain k = 3 . .
36 for shear stress and pressure,respectively. Since k is the amplitude of the leadingcorrection-to-scaling term, which is supposed to be asmall term, k must be O (1). This dramatically dependson the window of ˙ γ . If this window is far from the criticalregion, then the next terms in the correction-to-scalingmust be considered. Moreover, for such cases where thewindow of ˙ γ is far from the critical region and only theleading correction-to-scaling is considered, the obtainedvalue of k becomes too small or too large. Here, wesee that we arrive at conclusive values of k ∼ O (1)for both σ and p . This consistency check is crucial forthe analysis and must be carried out to examine thepre-assumptions for the correction-to-scaling terms.As a final note, the exponent y describes how the yieldstress σ y scales with distance to jamming δφ . This expo- -6 -5 -4 a m = d l n σ xy / d l n γ ˙ shear rate -6 -5 -4 b m = d l n P / d l n γ ˙ shear rate
50 0.843100 0.843150 0.843200 0.843250 0.843 50 0.844100 0.844150 0.844200 0.844250 0.844 300 0.8433 50 0.84335 100 0.84335 150 0.84335 200 0.84335 250 0.84335300 0.84335, 7300 0.8434300 0.84345
FIG. 3.
Successive slope of refined measurements.
Successive slope versus shear rate for ( a ) shear stress and ( b ) pressure.The critical density is found φ J = 0 . L ≤ L ≥ L = 300, the system consists of N = 65397 particles. The successiveslopes at slightly above and below φ J bend in opposite directions. A stronger correction-to-scaling is found for the shear stress σ . nent can be measured by simulations of pure isotropiccompression, and no-shearing is required. It is wellknown that y ≃ V. DISCUSSION AND CONCLUSION
Soft spheres flow freely in a dilute regime and becomeamorphous solid in a dense regime. This accounts fora large range of phenomena such as jamming and glasstransition. Determinations of both the transition densityand exponents describing scaling near the transitionpoint are subjects of intense research. However, becauseof a lack of a general framework, no consensus has yetemerged. Here, we close this debate by presenting aframework to precisely compute the exponents of therigidity transition in soft spheres based on an accuratedetermination of the transition density. Furthermore,we demonstrate that even though the transition densitycan be uniquely determined, there is a spectrum ofdifferent numerical values for critical exponents. Thanksto isotropic asymptotic scaling of the components of thestress tensor, we introduce a matching rule that selectscritical exponents and rules out other possibilities. Thematching rule considers the intersection of contour linesof exponents of pressure and shear stress. This allowedus to unambiguously determine the asymptotic criticalexponent of the shear stress and pressure. Havingdetermined the asymptotic exponent q , we use finite-sizescaling to determine the asymptotic value of exponentof the leading correction-to-scaling term ω/z . We demonstrate that ω/z for both shear stress and pressureconverges to the same value within the numericaluncertainty at the limit of large system sizes. Twomean-field type calculations for the exponents of therigidity transition are proposed by Otsuki-Hayakawa [42]and DeGiuli et al. [17]. Our results for exponents arecloser to the predictions by the former.Noticeably, we recover inertial-Bagnold scaling p, σ ∝ ˙ γ at ˙ γ → particles, are considered. Ourresults indicate that such small system sizes are stronglyinfluenced by finite size effects.Our framework provides grounds for several immedi-ate investigations that will deepen our understanding of a q ω /z σ xy p b ω / z L -1 σ xy p FIG. 4.
Matching rule. ( a ) Contour lines of q and ω/z for p and σ . The matching rule selects critical exponents where thecontour lines of p and σ intersect. ( b ) We hold fixed q ( L = ∞ ) = 0 .
41 and obtain correction-to-scaling exponent ω/z via fittingthrough Eq. 7. The obtained value of ω/z ( L ) is plotted against L − . One can see that ω/z saturates for large system sizes to ω/z ( L = ∞ ) = 0 . q ω/zσ, p . . Exponents.
Numerical values of critical exponents. amorphous materials using rheology as the main tool:(i) Critical exponents of a phase transition can be in-fluenced by fluctuations and thus the dimension-ality of the system. However, these exponents donot significantly change above a critical dimension,known as the upper critical dimension d UC . Belowthis dimension, fluctuations are important. Above d UC , fluctuations are washed out and critical expo-nents are equal to the mean-field exponents. Theexact determination of d UC for the jamming tran-sition has been a challenge: the absence of a meanfield theory and the lack of a framework for theprecise measurement of critical exponents can beconsidered as the main reasons. Many authorshave suggested that d UC = 2 and that logarithmiccorrections-to-scaling are involved [28, 41, 43–45].The main reason for this is that critical exponentsappear to be the same for d = 2 and 3. Collapseof the data has been used widely in these studies tomeasure critical exponents. Our general approachcan be easily applied in accurately measuring criti-cal exponents in three dimensions. Then, a compar-ison of critical exponents at d = 2 and 3 can resolvethe controversy over the upper critical dimensionfor jamming. This will be a great step ahead in un- derstanding the nature of the jamming transition.(ii) Amorphous solids possess a complex free-energylandscape [46]. As one increases density, an amor-phous solid undergoes a sequence of transitions:glass transition, Gardner transition, and jamming.Annealing has been the essential method to investi-gate these transitions. Standard rheological tech-niques have been shown to be powerful tools toinvestigate complex properties of this energy land-scape of amorphous materials [47]. We expect thegeneralization of our approach to shed light on andhelp formulate a general formalism to investigateother types of transitions in amorphous materialsusing rheology.(iii) An investigation of periodically driven colloidal sus-pensions provided remarkable insights into the na-ture of rheological phase transitions. In a diluteregime, these systems undergo a non-equilibriumphase transition into an absorbing state whereparticles self-organize themselves to prevent colli-sions [48, 49]. In the dense regime, a yielding tran-sition, which describes the onset of plastic deforma-tion, has been shown to be a non-equilibrium phasetransition from reversibility in an elastic regime into irreversibility in a plastic regime [50]. Nonetheless,the nature of the transition has been disputed, in-cluding whether it is a first-order or second-ordertransition [51, 52], as well as whether the absorbingstate transition belongs to the universality class of(conserved) directed percolation [50]. In both thedilute and dense regimes, divergences of time andlength scales have been reported upon approachinga critical shearing strain. However, a precise de-termination of the universality classes of these non-equilibrium phase transitions have not been conclu-sive due to the lack of a comprehensive frameworkfor measuring critical exponents from rheologicalexperiments [48–50, 53, 54]. Our formalism mayshed light on resolving the dispute over the natureof rheological phase transitions in oscillatory shear-ing.Rheological phase transitions are fascinating noveltransitions, and the exploration of their characteristicsprovides new insights into the less-explored realm ofathermal non-equilibrium phase transitions [55, 56].Compared to other well-established equilibrium transi-tions, rheological phase transitions are in their infancy.We hope that our framework can aid in a better under-standing of their nature. Acknowledgments.
The authors thank the KoreaInstitute for Advanced Study for providing computingresources (KIAS Center for Advanced Computation- Linux cluster system) for this work, and especiallyconsultations from Hoyoung Kim. We appreciate en-lightening discussions with Abbas Ali Saberi, TakahiroHatano, Peter Olsson, and Hisao Hayakawa. Thiswork is supported in part by the NRF grant No.2017R1D1A1B06035497.
VI. APPENDIX A: SIMULATIONS
Numerical simulations.
We perform constant vol-ume molecular dynamics simulations of two-dimensionalfrictionless bidisperse disks. Interactions betweenparticles are modeled by a linear dashpot model. Twoparticles i and j of radii R ai and R bj (where a, b = 0 and1 stand for two different radii of bidisperse particles) atpositions r i and r j interact when ξ ij = R ai + R bj − r ij > ξ ij is called the mutual compression of particles i and j , r ij = | r i − r j | . The particles interact via a lineardashpot model, F ij = Y ξ ij + γ dξ ij dt , where Y and γ aredenoted as elastic and dissipative constants, respectively.Throughout the study, we adopt unitary scale Y = 1and γ = 1, respectively.To prevent crystallization, we use a 1 : 1 binarymixture of particles where the ratio of the radii of largeand small particles is set to R /R = 1 .
4. The diameter of small particles is chosen as the unit of the length2 R = 1, and the mass of each particle is equal to itsarea, m a = π [ R a ] .Lees-Edwards boundary conditions are applied alongthe y -direction. They create a uniform overall shearrate, ˙ γ . We use LAMMPS for our simulations. Thanksto the developer team of LAMMPS, we were providedwith a new version of LAMMPS that prevents artificialattractive forces arising from the dashpot model. Theversion can be accessed via the mailing list of LAMMPS.We used several system sizes, the smallest L = 50 andthe largest L = 300. We change the packing fraction bychanging the number of particles N via: N = 8 π L + 1 . φ, (10)where 1 and 1 . γ in the range10 − and 10 − , the total strain is γ = 30 L and theintegration time step is dt = 0 .
1. For the next smallerdecade, the integration time step is dt = 0 . VII. APPENDIX B: SCALING ANSATZ
Here, we explain a formalism for deriving the scalingansatz for a rigidity transition. The formalism in prin-ciple can be applied to any transition that is accompa-nied by a diverging length scale ξ . Upon approachingthe dense regime, the motion of particles becomes co-ordinated. This signals the growing length scale, whichdiverges at the critical density φ J . This divergence isdescribed by exponent ν via: ξ ∝ δφ − ν . (11)In the proximity of a critical point, the only fundamentallength scale b is the correlation length scale, b = ξ . Eq. 11can be cast into a dimensionless number asΠ φ = δφb /ν . (12)The critical point is at δφ = 0 and ˙ γ →
0, therefore at δφ = 0, the correlation length diverges upon decreasingthe shear rate: ξ ∝ ˙ γ − /z , (13)where z is the dynamic exponent. This equation can besimilarly cast into another dimensionless number viaΠ ˙ γ = ˙ γb z . (14)Now, any physical quantity such the shear stress σ alsoscales with the distance from jamming σ ∝ δφ y at ˙ γ → σ ∝ b − y/ν , (15)which provides the dimensionless number for this quan-tity Π σ = σb y/ν . (16)Since σ depends on both δφ and ˙ γ Π σ = F (Π δφ , Π ˙ γ ) , (17)which results in σb y/ν = F (cid:16) δφb /ν , ˙ γb z (cid:17) , (18)This is the dimensionless equation of state.In the renormalization group method, the domainover the correlated particles are rescaled. After renor-malization, the system becomes smaller by a factor of b and therefore ξ ′ = ξ/b . As a result of this, the systemmoves away from the critical point by renormalization.In this process, all observables and control parametersscale with distance from critical point b . Equation 18describes all such scaling behaviors. Two approaches,the intermediate asymptotic approach described bydimensionless numbers and the renormalization group,arrive at similar results [57].If we choose the length scale b such that ˙ γb z = 1, then σ = ˙ γ y/zν F (cid:18) δφ ˙ γ /zν (cid:19) , (19)which is the leading scaling term. At δφ = 0, σ ∝ ˙ γ q ,thus q/y = 1 /zν . This equation describes σ infinitesi-mally close to the critical point at δφ = 0 and ˙ γ = 0. The jamming point is characterized by two principaldirections given by δφ and ˙ γ . Each direction is accompa-nied by a principal exponent: y and q . Near the criticalpoint only these relevant quantities affects the dynamics.However, off the critical point, some irrelevant parame-ters, w , may affect the dynamics. Since this quantity isirrelevant, one cannot bring the system into the criticalpoint by varying such a quantity. This means that thecorrelation length does not diverge if w →
0. However,it may retain a scaling form near the critical region ξ ∝ w /ω , (20)which results to Π w = wb − ω . (21)Inserting this dimensionless number into Eq. 18 resultsin σb y/ν = F (cid:16) δφb /ν , ˙ γb z , wb − ω (cid:17) . (22)With ˙ γb z = 1, we arrive at σ = ˙ γ y/zν F (cid:18) δφ ˙ γ /zν , w ˙ γ ω/z (cid:19) , (23)A Taylor expansion of this equation to the first ordergives: σ = ˙ γ y/zν (cid:20) F (0)1 (cid:18) δφ ˙ γ /zν (cid:19) + ˙ γ ω/z F (1)1 (cid:18) δφ ˙ γ /zν (cid:19)(cid:21) , (24)This equation describes the leading correction-to-scalingterm. At δφ = 0 σ = ˙ γ y/zν h c + c ˙ γ ω/z i . (25)Eq. 25 can be used for scalings of the flow curve at φ J . [1] D. Bonn, M. M. Denn, L. Berthier, T. Divoux, andS. Manneville, “Yield stress materials in soft condensedmatter,” Rev. Mod. Phys. , 035005 (2017).[2] G. Petekidis, D. Vlassopoulos, and P. N. Pusey, “Yield-ing and flow of sheared colloidal glasses,” J. Phys.: Con-dens. Matter , S3955 (2004).[3] A. J. Liu and S. R. Nagel, “Nonlinear dynamics: Jam-ming is not just cool any more,” Nature , 21 (1998).[4] J. Paredes, M. A. J. Michels, and D. Bonn, “Rheologyacross the zero-temperature jamming transition,” Phys.Rev. Lett. , 015701 (2013).[5] V. Trappe, V. Prasad, L. Cipelletti, P. N. Segre, andD. A. Weitz, “jamming phase diagram for attractive par-ticles,” Nature , 772 (2001).[6] S. H. E. Rahbari, M. Khadem-Maaref, and S. K.A. Seyed Yaghoubi, “Universal features of the jammingphase diagram of wet granular materials,” Phys. Rev. E , 042203 (2013).[7] H. A. Barnes and K. Walters, “The yield stress myth?”Rheologica acta , 323–326 (1985).[8] P. C. F. Moeller, A. Fall, and D. Bonn, “Origin of ap-parent viscosity in yield stress fluids below yielding,” Eu-rophys. Lett. , 38004 (2009).[9] P. Olsson and S. Teitel, “Critical scaling of shearing rhe-ology at the jamming transition of soft-core frictionlessdisks,” Phys. Rev. E , 030302 (2011).[10] E. Lerner, G. D¨uring, and M. Wyart, “A unified frame-work for non-Brownian suspension flows and soft amor-phous solids,” Proc. Nat. Acad. Sc. , 4798 (2012).[11] T. Kawasaki, D. Coslovich, A. Ikeda, and L. Berthier,“Diverging viscosity and soft granular rheology in non-Brownian suspensions,” Phys. Rev. E , 012203 (2015).[12] P. Schall and M. van Hecke, “Shear bands in matter withgranularity,” Annu. Rev. Fluid Mech. (2010). [13] T. Hatano, “Scaling properties of granular rheology nearthe jamming transition,” J. Phys. Soc. Jpn. , 12 (2008).[14] M. Otsuki and H. Hayakawa, “Universal scaling for thejamming transition,” Prog. of Theor. Phys. , 647–655(2009).[15] T. Hatano, “Critical scaling of granular rheology,” Progr.Theor. Exp. Phys. , 143 (2010).[16] M. Otsuki and H. Hayakawa, “Rheology of sheared gran-ular particles near jamming transition,” Prog. of Theor.Phys. Supp. , 129–138 (2012).[17] E. DeGiuli, G. D¨uring, E. Lerner, and M. Wyart, “Uni-fied theory of inertial granular flows and non-Browniansuspensions,” Phys. Rev. E , 062206 (2015).[18] D. Vagberg, P. Olsson, and S. Teitel, “Critical scalingof Bagnold rheology at the jamming transition of fric-tionless two-dimensional disks,” Phys. Rev. E , 052902(2016).[19] D. Vagberg, P. Olsson, and S. Teitel, “Dissipation andrheology of sheared soft-core frictionless disks below jam-ming,” Phys. Rev. Lett. , 208303 (2014).[20] P. Olsson and S. Teitel, “Herschel-bulkley shearing rhe-ology near the athermal jamming transition,” Phys. Rev.Lett. , 108001 (2012).[21] D. Vagberg, P. Olsson, and S. Teitel, “Universality ofjamming criticality in overdamped shear-driven friction-less disks,” Phys. Rev. Lett. , 148002 (2014).[22] F. Radjai and S. Rou, “Turbulentlike fluctuations in qua-sistatic flow of granular media,” Phys. Rev. Lett. ,064302 (2002).[23] P. Olsson, “Relaxation times and rheology in dense ather-mal suspensions,” Phys. Rev. E , 062209 (2015).[24] K. N. Nordstrom, E. Verneuil, P. E. Arratia, A. Basu,Z. Zhang, A. G. Yodh, J. P. Gollub, and D. J. Durian,“Microfluidic rheology of soft colloids above and belowjamming,” Phys. Rev. Lett. , 175701 (2010).[25] P. Olsson and S. Teitel, “Critical scaling of shear viscosityat the jamming transition,” Phys. Rev. Lett. , 178001(2007).[26] M. Kardar, Statistical physics of fields (Cambridge Uni-versity Press, 2007).[27] H. Hayakawa and M. Otsuki, “Critical behaviors ofsheared frictionless granular materials near the jammingtransition,” Phys. Rev. E , 032117 (2013).[28] C. P. Goodrich, A. J. Liu, and J. P. Sethna, “Scalingansatz for the jamming transition,” Proc. Nat. Acad. Sc. , 9745 (2016).[29] C. E. Maloney and A. Lemaˆıtre, “Amorphous systems inathermal, quasistatic shear,” Phys. Rev. E , 016118(2006).[30] E. Bouchbinder, J. S. Langer, and I. Procaccia, “Ather-mal shear-transformation-zone theory of amorphous plas-tic deformation. i. basic principles,” Phys. Rev. E ,036107 (2007).[31] A. Lemaitre and C. Caroli, “Rate-dependent avalanchesize in athermally sheared amorphous solids,” Phys. Rev.Lett. , 065501 (2009).[32] H. G. E. Hentschel, S. Karmakar, E. Lerner, and I. Pro-caccia, “Size of plastic events in strained amorphoussolids at finite temperatures,” Phys. Rev. Lett. ,025501 (2010).[33] T. Hatano, C. Narteau, and P. Shebalin, “Com-mon dependence on stress for the statistics of granularavalanches and earthquakes,” Sci. Rep. (2015).[34] A. Nicolas, E. E. Ferrero, K. Martens, and J. L. Barrat, “Deformation and flow of amorphous solids: a reviewof mesoscale elastoplastic models,” arXiv:1708.09194(2017).[35] L. D. Landau and E. M. Lifshitz, “Theory of elasticity,vol. 7,” Course of Theoretical Physics , 109 (1986).[36] P. E. Peyneau and J. N. Roux, “Solidlike behavior andanisotropy in rigid frictionless bead assemblies,” Phys.Rev. E , 041307 (2008).[37] M. Baity-Jesi, C. P. Goodrich, A. J. Liu, S. R. Nagel, andJ. P. Sethna, “Emergent so3 symmetry of the frictionlessshear jamming transition,” J. Stat. Phys. , 735–748(2017).[38] D. Vagberg, P. Olsson, and S. Teitel, “Effect of collisionalelasticity on the bagnold rheology of sheared frictionlesstwo-dimensional disks,” Phys. Rev. E , 012902 (2017).[39] K. Suzuki and H. Hayakawa, “Theory for the rheology ofdense non-brownian suspensions: divergence of viscosi-ties and µ - j rheology,” arXiv:1711.08855 (2017).[40] C. Heussinger and J. L. Barrat, “Jamming transition asprobed by quasistatic shear flow,” Phys. Rev. Lett. ,218303 (2009).[41] C.S. O’Hern, L.E. Silbert, A.J. Liu, and S.R. Nagel,“Jamming at zero temperature and zero applied stress:The epitome of disorder,” Phys. Rev. E , 011306(2003).[42] M. Otsuki and H. Hayakawa, “Critical behaviors ofsheared frictionless granular materials near the jammingtransition,” Phys. Rev. E , 011308 (2009).[43] M. Wyart, L. E. Silbert, S. R. Nagel, and T. A. Wit-ten, “Effects of compression on the vibrational modesof marginally jammed solids,” Phys. Rev. E , 051306(2005).[44] C. P. Goodrich, A. J. Liu, and S. R. Nagel, “Finite-sizescaling at the jamming transition,” Phys. Rev. Lett. ,095704 (2012).[45] C. P. Goodrich, S. Dagois-Bohy, B. P. Tighe, M. vanHecke, A. J. Liu, and S. R. Nagel, “Jamming in finitesystems: Stability, anisotropy, fluctuations, and scaling,”Phys. Rev. E , 022138 (2014).[46] P. Charbonneau, J. Kurchan, G. Parisi, P. Urbani, andF. Zamponi, “Fractal free energy landscapes in structuralglasses,” Nat. Commun. , 3725 (2014).[47] Y. Jin and H. Yoshino, “Exploring the complex free-energy landscape of the simplest glass by rheology,” Nat.Commun. , 14935 (2017).[48] D. J. Pine, J. P. Gollub, J. F. Brady, and A. M. Leshan-sky, “Chaos and threshold for irreversibility in shearedsuspensions,” Nature , 997 (2005).[49] L. Corte, P. M. Chaikin, J. P. Gollub, and D. J. Pine,“Random organization in periodically driven systems,”Nature Phys. , 420 (2008).[50] K. H. Nagamanasa, S. Gokhale, A. K. Sood, andR. Ganapathy, “Experimental signatures of a nonequi-librium phase transition governing the yielding of a softglass,” Phys. Rev. E , 062308 (2014).[51] I. Regev, J. Weber, C. Reichhardt, K. A Dahmen, andT. Lookman, “Reversibility and criticality in amorphoussolids,” Nat. Commun. , 8805 (2015).[52] P. Leishangthem, A. D. S. Parmar, and S. Sastry, “Theyielding transition in amorphous solids under oscillatoryshear deformation,” Nat. Commun. , 14653 (2017).[53] E. Tjhung and L. Berthier, “Hyperuniform density fluc-tuations and diverging dynamic correlations in periodi-cally driven colloidal suspensions,” Phys. Rev. Lett. , , 3474 (2014).[55] Haye Hinrichsen, “Non-equilibrium critical phenomenaand phase transitions into absorbing states,” Adv. Phys. , 815 (2000).[56] Z. R´acz, “Nonequilibrium phase transitions,” arXivpreprint cond-mat/0210435 (2002).[57] N. Goldenfeld, O. Martin, and Y. Oono, “Intermediateasymptotics and renormalization group theory,” J. Sci.Comput.4