A direct link between active matter and sheared granular systems
Peter K. Morse, Sudeshna Roy, Elisabeth Agoritsas, Ethan Stanifer, Eric I. Corwin, M. Lisa Manning
AA direct link between active matter and sheared granular systems
Peter K. Morse, ∗ Sudeshna Roy, ∗ Elisabeth Agoritsas, ∗ Ethan Stanifer, Eric I. Corwin, and M. Lisa Manning Department of Chemistry, Duke University, Durham, North Carolina 27710, USA Department of Physics, Syracuse University, Syracuse, New York 13244, USA Institute of Physics, Ecole Polytechnique F´ed´erale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland Department of Physics and Materials Science Institute, University of Oregon, Oregon 97403, USA (Dated: September 17, 2020)The similarity in mechanical properties of dense active matter and sheared amorphous solidshas been noted in recent years without a rigorous examination of the underlying mechanism. Wedevelop a mean-field model that predicts that their critical behavior should be equivalent in infinitedimensions, up to a rescaling factor that depends on the correlation length of the applied field.We test these predictions in 2d using a new numerical protocol, termed ‘athermal quasi-staticrandom displacement’, and find that these mean-field predictions are surprisingly accurate in lowdimensions. We identify a general class of perturbations that smoothly interpolate between theuncorrelated localized forces that occur in the high-persistence limit of dense active matter, andsystem-spanning correlated displacements that occur under applied shear. These results suggest auniversal framework for predicting flow, deformation, and failure in active and sheared disorderedmaterials.
SIGNIFICANCE STATEMENT
There is not yet a robust theoretical frameworkpredicting the dynamics of dense active matter,where energy is injected at the scale of constituentparticles. Previous work has shown that some fea-tures of dense active matter are similar to thosein dense disordered materials that are shearedglobally from the boundaries. Using analytic andcomputational tools, we show that there is a di-rect correspondence between active matter andapplied shear strain, which can in turn be usedto help predict the behavior of dense active mat-ter.The statistical physics of active matter – where en-ergy is injected at the smallest scale, that of the particlesthemselves – is highly nontrivial, exhibiting new featuressuch as giant number fluctuations and motility-inducedphase separation [1, 2]. While comprehensive theorieshave been developed for many of these phenomena at lowand intermediate densities [2, 3], the behavior of highlydense, glassy active matter remains more mysterious. Re-cent work by Henkes and collaborators [4, 5] highlightsthe important role of the potential energy landscape inconstraining and dictating the behavior of dense activematter, which is in some ways similar to the situationin glasses excited by thermal fluctuations. Nevertheless,work by Berthier and collaborators emphasizes impor-tant differences between the dynamics of thermal andactive glasses [6, 7] within the glassy potential energylandscape. Therefore, the large body of work on ther-mally excited glasses can not be transferred immediately ∗ These authors contributed equally to this work to active glasses, and so a predictive theory for the dy-namics of dense active matter remains elusive.Meanwhile, the dynamics of athermal sheared disor-dered materials, where energy is injected at the largestscale, globally from the boundaries, has been the sub-ject of intense study for decades. A recent breakthroughallows an exact analytic solution for the behavior ofslowly sheared systems in infinite dimensions, where in-teractions are exactly mean-field [8–12]. These resultsqualitatively explain many features in sheared 2- and 3-dimensional glassy solids. Perhaps more interestingly,new work suggests that the dynamical mean-field equa-tions in infinite dimensions have the same structure re-gardless if the driving forces are generated by global shearor active forces on each particle [13, 14], as all such forc-ing can be represented by memory kernels with the samefunctional form.There is also evidence of similarities between shearedand active glassy systems in 2- and 3-dimensional simula-tions; recent studies have noted that in granular systemsthe two forcing mechanisms yield viscosities [15] and largedensity fluctuations [4, 16] that are similar.What is missing in the low-dimensional scenarios is aunifying picture as developed in infinite dimensions; todevelop such a picture, it is necessary to first examinehow and where discrepancies between shear and randomforces appear. For example, Liao and Xu [15] noted thatself-propelled particles driven by constant forces with thesame magnitude in random directions will have the samedynamic viscosity as their sheared counterparts [17–19],albeit with different critical exponents. Moreover, thevalues of the exponents can be changed by altering fea-tures of the forces on the self-propelled particles. There-fore, one wonders whether there may be a family offorcing fields, including shear and different types of self-propulsion, where all the resulting dynamics could be un-derstood and predicted as part of a universal descriptionof failure in jammed solids. a r X i v : . [ c ond - m a t . s o f t ] S e p One hint about how such a framework might be con-structed comes from the density of states that describesthe spectrum of vibrational modes about a mechanicallystable state in the potential energy landscape. Morespecifically, in low dimensions, it has been shown thatthe linear response of particles to either random forcesin the limit of low rotational noise or long persistencelength [4, 5, 20], or to shear [21] is dominated by thelowest eigenmode. Very close to an instability, this low-est eigenmode specifies the direction in the energy land-scape with the lowest energy barrier [22], and highlightsthe direction in which particles must move to leave onemechanically stable state and find another [23, 24].Taken together, these previous results suggest that in2d and 3d materials there is a direct connection betweenhow a disordered system traverses the energy landscapeunder shear and under random forces in the limit ofzero rotational noise. Here we develop an exact infinite-dimensional mean-field theory prediction for the mechan-ical response of materials under shear and active forces.We explicitly test this prediction by analyzing numeri-cal simulations of soft spheres in two dimensions, andcomparing dynamics under athermal quasi-static shear(AQS) [24] and a new constrained dynamics we termathermal quasi-static random displacements (AQRD).The latter is equivalent to self-propelled forcing in thelimit of zero rotational noise, where the forcing is slowerthan any other relaxation process inside the material.We show that under these two different types of ather-mal driving, scaling relations (including exponents) de-scribing the avalanche statistics and the sampling of sad-dle points are identical and consistent with mean-fieldpredictions, although the coefficients of these scaling re-lations differ.We conjecture that differences in those coefficients, in-cluding the shear modulus, are governed by the corre-lation lengthscale associated with the imposed displace-ment field; in shear this length is the size of the box,while for completely random fields it is the size of indi-vidual particles. In addition, the mean-field calculationpredicts that the coefficients should scale with the fluc-tuations in strain between nearby particles, which alsochange systematically with the correlation length of theimposed displacement field.Therefore, we systematically vary this correlationlength in our simulations and find that the coefficients ex-hibit a systematic power-law scaling that matches mean-field predictions. We also study the effect of materialpreparation on these results, demonstrating that shearand random displacement fields are similar even in ultra-stable glasses. This demonstrates that shear can be con-sidered as a highly-correlated special case of more generalrandom displacements.
I. METHODS
Two ways of traversing the energy landscape—
When constructing the energy landscape of allowed configura-tions, there are two types of variables that play a pri-ori different roles: state variables are explicitly speci-fied by the experimental or simulation protocol, while reaction coordinates are free to vary under constraintsimposed by state variables. For instance, a standardinfinite temperature quench [25] considers shear-strainto be a state variable during preparation, while shear-stabilization methods [26] treat strain as a reaction co-ordinate during preparation, regardless of how the strainvariable is used afterwards. Therefore, the use of strainor the box degrees of freedom as state variables is merelyan artifact of the way in which experiments or simula-tions are performed. Moreover, during an athermal qua-sistatic perturbation, we adjust a state variable and thenre-minimize the system by allowing all reaction coordi-nates to find their nearest local energy minima.An applied shear strain, illustrated by the red arrowsin Fig. 1a, perturbs the system in its the potential energylandscape. One way to represent this perturbation is toview the landscape as a function of the
N d reaction co-ordinates (particle positions), so that adjusting the statevariable (shear strain) contorts the landscape in that
N d -dimensional space [23, 24, 26]. As a system is shearedtowards a saddle point, a nearby energy barrier is low-ered until the system reaches the saddle point and movesdownhill towards a new minimum. It is equivalent to de-scribe this process instead as moving along an
N d + N d -dimensional representation, whereasthe saddles parallel to the strain correspond to the shearmodulus changing sign, which does not correspond to aninstability in a strain-controlled measurement [27].A second type of possible perturbation is a randomdisplacement field, where we choose a random directionin configuration space ∣ c ⟩ and promote it to a state vari-able. An example field ∣ c ⟩ is illustrated by the red arrowsin Fig. 1d. Thus, after perturbing along ∣ c ⟩ , the systemis free to relax along all directions perpendicular to ∣ c ⟩ ,but motion along ∣ c ⟩ is restricted via constrained min-imization to the other N d − ∣ c ⟩ , and we ask whether the distribution of saddlesand their corresponding stress drops follow the same dis-tribution as those encountered under shear strain. Numerical Model Description—
We simulate N Hertzian spheres in d = N is thenumber of particles. Except where specified when usingultrastable glasses, our systems are a 50-50 mixture ofbidisperse disks with diameter ratio 1:1.4 to avoid crys-talization. For the pressure sweep data, we prepare oursystems at a target pressure by performing a standardinfinite temperature quench [25], followed by FIRE min-imization [28] at a packing fraction such that we stay Reaction Coordinate(1 of Nd) U S t r a i n b -4 d U Reaction Coordinate(1 of Nd - 1) eac f -3 FIG. 1.
Two methods of traversing the energy land-scape: AQS and AQRD . a) Forces applied to particlesin an AQS ensemble. b) Potential energy landscape splittingout the Nd + Nd (one of which is thereaction coordinate shown) and strain. c) Stress-strain curveshowing that stress drop occurs when a saddle point in thereaction coordinate is reached by traversing along the straincoordinate. d) Forces applied to particles in a sample AQRDensemble. e) Potential energy landscape splitting out the Nd degrees of freedom (for fixed box shape) into Nd − f ) Random-stressvs. random-strain curve showing that random-stress drop oc-curs when a saddle point in the reaction coordinate is reachedby traversing along the ∣ c ⟩ coordinate. above the target pressure, followed by a careful decom-pression [27, 29]. For the correlation length sweep, weprepare our systems at a pressure of p = . ± . φ = .
94 [25]. In each case, we use the Hertziancontact potential U = / ∑ ij Θ ( ε ij ) ε / ij (1)where Θ is the Heaviside function, ε ij = − r ij /( ρ i + ρ j ) is the dimensionless overlap, ρ i is the radius of particle i , and r ij is the distance between particles i and j . Alllength scales are reported in natural units of the averageparticle diameter. Athermal Quasi-static Shear—
Under the now-standard method of Athermal Quasi-static Shear(AQS) [24], our system of particles is subject to Lees-Edwards boundary conditions where the periodic replicas in the y -direction are shifted by an amount γL y in the x -direction, and γ is the shear strain. After each small stepin the applied strain (∆ γ = − ), a FIRE minimizationalgorithm [28] is used to minimize the energy subject tothe constraint that the box shape is held fixed. There-fore, AQS is equivalent to dynamics in the limit of zerostrain rate – where the material is sheared more slowlythan any process or relaxation rate inside the material.To facilitate comparison with the AQRD protocol de-scribed in the next section we emphasize that, in linearresponse and neglecting the effect of particle-particle in-teractions, shearing the boundary a distance γL y alongthe x -direction is equivalent to displacing particles inthe x -direction with a magnitude determined via theheight of the system as given by u αi = γδ αx ( y i − L y / ) .Here y i is the y -coordinate of particle i , L y is thelength of the box in the y -direction, and δ is the Kro-necker delta function of x and dimensional index α [24].An example of such a displacement field is shown inFig. 1a. The overall magnitude of this displacementvector field generated by an applied strain γ is thengiven by ∣ u ( γ )∣ = γ [ ∑ i ( y i − L y / ) ] / . If we assume auniform distribution of y -coordinate values, as one ex-pects in an amorphous sample, the average magnitude is ∣ u ( γ )∣ ≈ γL y √ N /
12. Therefore, an applied shear strainof γ is equivalent to moving a distance γL y √ N /
12 alonga normalized vector field, independent of dimension.
Athermal Quasi-static Random Displacements—
Sim-ilar to AQS, the system is initialized into a mechanicallystable state at the bottom of a potential energy well withenergy U and N d -dimensional position vector ∣ x min ⟩ .The system is then displaced along an N d -dimensionalunitless vector ∣ c ⟩ with elements c i and ⟨ c ∣ c ⟩ =
1. We ex-plore different methods for choosing ∣ c ⟩ described below.First, we define the random strain ˜ γ as a displacement ˜ u along the vector ∣ c ⟩ scaled as˜ γ = ˜ uL y √ N . (2)This definition ensures that strains ( γ ) in AQS can bedirectly compared to random-force strains (˜ γ ) in AQRD,where both are unitless.Starting from positions ∣ x min ⟩ and displacing by anamount ˜ u , new positions are then ∣ x ⟩ = ∣ x min ⟩ + ˜ u ∣ c ⟩ , butthey are not in a local energy minimum with respect tothe reaction coordinates. Therefore, we must evolve thesystem using a constrained minimization that imposesan external force ∣ F ext ⟩ = − λ ∣ c ⟩ , where λ is the Lagrangemultiplier, which prevents any motion along ∣ c ⟩ .We calculate how such displacements induce changesto the internal stress of the system, in direct analogy tostress-strain curves for AQS. The stress induced by thefield ˜ γ is given by˜ σ = dUd ˜ γ = N ∑ i = ( ∂U∂x i dx i d ˜ γ + ∂U∂x ⊥ i dx ⊥ i d ˜ γ ) . (3)where we have split the particle motion x i into com-ponents which are parallel or perpendicular to c i as x i and x ⊥ i respectively. By definition, ∂U∂x ⊥ i = − F ⊥ i andsince we minimize force with respect to the particle po-sition, F ⊥ i =
0. Thus, the total residual force F i on eachparticle i is parallel to c i . Furthermore, we note that dx i d ˜ γ = c i L y √ N , resulting in the definition of the randomstress ˜ σ = − N ∑ i = F i ⋅ c i L y √ N = − ⟨ F ∣ c ⟩ L y √ N . (4)This is a generalization of the derivation for shear stressin AQS developed Maloney and Lemaitre [24]. Through-out this manuscript, we use variables with a tilde to de-note observables that are the AQRD equivalent to AQScounterparts.In practice, we evolve the system by taking steps of10 − in the random strain ˜ γ , and after each step we useFIRE minimization [28] to find the constrained local min-imum. Thus, instead of applying forces in the FIRE-calculated gradient direction ∣ F ⟩ , we apply them along ∣ F ⟩ − ⟨ c ∣ F ⟩ ∣ c ⟩ . We impose a stopping condition when ev-ery component of the total excess force on every particleis less than a cutoff value of 10 − , set to ensure particlepositions to double precision. By construction, there isno drift velocity in the system.While this version of AQRD applies displacements ina direct analogy to a strain-controlled measurement, wedescribe a stress-controlled version of random forcing inthe supplement. We show that the two ensembles givethe same results in linear response. Thus, while AQRDapplies a random strain, this is equivalent to applying arandom body force on every particle, at least until thematerial yields.We generate the fields ∣ c ⟩ for AQRD using two differentmethods ways: one based on random Gaussian fields andanother based on plane waves. The Gaussian randomfields, which are spatially correlated over a characteris-tic length scale ξ , are generated using a standard Fouriertransform method that respects the periodic boundaryconditions. A detailed description is given in the supple-ment. Fig. 2a-c illustrates the random vector ∣ c ⟩ gener-ated from the correlated Gaussian random field for dif-ferent correlation length ξ =
1, 2 . .
25, respectively.To test whether features we observe are dependent onlyon the correlation length, or whether other features ofthe field structure are important, we also generate plane-wave-like fields where the x − components of the vectorsare a sine function of the y − coordinate of the particle po-sitions, and the y − components of the vectors vanish. Forsuch fields, we define the correlation length scale to behalf the chosen wavelength. Fig. 2e-g illustrates the ran-dom vector ∣ c ⟩ generated in wave-like pattern for differentcorrelation lengths ξ = .
5, 6 .
25 and 25 respectively. Thecorresponding displacement field associated with shearunder Lees-Edwards boundary conditions is equivalentto a plane wave with a wavelength 2 L y , which is clear from Fig. 1a. II. RESULTS
Mean-field results—
The limit of infinite dimension pro-vides an exact benchmark to investigate properties ofstructural glasses [12, 30], and has been successfullyused, for instance, to study quasistatic shear or compres-sion [8, 11]. In this framework, we can show that AQSand AQRD are strictly equivalent upon a simple rescal-ing of the accumulated strain, with a dependence on thecorrelation length ξ . The full derivation is provided inRef. [31].In order to implement a local strain vector ∣ c ⟩ ∈ R Nd as in AQRD, we assign to each particle a random localstrain c i drawn from a Gaussian distribution with zeromean defined by: c i = , c i c j = Ξ f ξ (∣ r ij ( )∣)/ d , with f ξ ( x ) = e − x /( ξ ) /√ πξ , (5)where the overline denotes the statistical average overthe quenched random strain field, Ξ is a tunable ampli-tude which has the units of a length (so that the strainsremain unitless), and r ij ( t ) is the distance between par-ticles i and j at time t (and we focus here on the initialconfiguration). For simplicity we have assumed the fluc-tuations in the field can be described by a normalizedGaussian function with a finite correlation length ξ > d so thatthe fluctuations in c scale with dimension in the sameway as fluctuations in the local strain field in AQS.In the infinite-dimensional limit, the complex many-body dynamics of pairwise interacting particles becomesexactly mean-field and can be reduced to an effectivescalar stochastic process for the fluctuating gap betweenparticle pairs, h ij ( t ) = d (∣ r ij ( t )∣/ (cid:96) − ) where (cid:96) is the typ-ical distance between particles and h ij ∼ O( ) [13, 14, 32].The dynamics are governed by the distribution of the rel-ative strains c ij ≡ c i − c j , which are uncorrelated in thelimit d → ∞ for distinct pairs of particles (consistent withthe mean-field assumption). The variance of a given pair,however, still encodes the spatial correlation of individuallocal strains, through the following quantity: F ( Ξ , (cid:96), ξ ) = d(cid:96) c ij = (cid:96) Ξ [ f ξ ( ) − f ξ ( (cid:96) )] , (6)which can be straightforwardly computed for a givenchoice of f ξ , or directly measured in numerical simula-tions. By adapting the derivation of the mean-field de-scription for shear presented in Ref. [13], we find thatAQS and AQRD are equivalent in infinite dimension, pro-vided that we rescale the accumulated strain by a factor √ F / (cid:96) , so that it is directly controlled by the variance ofrelative strains c ij .For the quasistatic stress-strain curves and the elas-tic modulus, we specifically predict that the random FIG. 2.
Effect of random field correlation length ξ on the mechanical response. a-c Snapshots of Gaussian correlatedfields (GCF) with correlation lengths a) ξ = b) ξ = . c) ξ = . d) Example random-stress vs. random-strain curvesfor random fields with different correlation lengths. e-g
Snapshots of wave-like correlated fields (WCF) with wave lengths e) ξ = . f ) ξ = . g) ξ = h) Example random-stress vs. random-strain curves for wave-like fields with differentcorrelation lengths. strain ˜ γ can be written in terms of the AQS shear strain γ , and therefore the random-displacement stress ˜ σ andthe random-displacement modulus ˜ µ can also be easilyscaled: γ MF ≡ ˜ γ MF √ F (cid:96) ⇒ ⎧⎪⎪⎨⎪⎪⎩ σ MF = (cid:96) √ F ˜ σ MF ,µ MF = (cid:96) F ˜ µ MF , (7)where the MF subscripts emphasize that this is a mean-field prediction, whose validity should be tested in lowerdimensions.We emphasize that the infinite-dimensional calculationpredicts that F / (cid:96) is thus the key quantity to make theAQRD random stress-random strain curves (and othersuch mean-field observables) collapse onto their AQScounterparts, as the dynamics are completely governedby the variance of relative local strains on particle pairs.Under our assumption that f ξ ( x ) is a normalizedGaussian function as in Eq. (5), we can straightforwardlycompute F from Eq. (6). By Taylor-expanding F in thelimits (cid:96) / ξ ≪ (cid:96) / ξ ≫
1, and keeping only the leadingterms, we predict a crossover of the elastic modulus ξ -dependence depending on the ratio (cid:96) / ξ , with F ∼ / ξ at (cid:96) / ξ ≫ F ∼ / ξ at (cid:96) / ξ ≪ ξ is of the order of the system size for shear. In bothcases, this implies that the elastic modulus decreases withincreasing ξ , as we will demonstrate numerically below,and matches with physical intuition: it is easier to de- form a glass with more coordinated local strains, i.e. witha larger correlation length. In particular, Eq. 7 statesthat in mean-field, AQS is a special case of AQRD, with F / (cid:96) = Numerical results for random stress vs. random strain–
We next test the mean-field prediction in numericalsimulations in 2D. Our first observation is that AQS andAQRD give rise to qualitatively similar stress-strain andrandom stress-random strain curves, as highlighted inFigs. 1c and 1f. Elastic branches–where the stress riseslinearly with the strain–are punctuated by points wherethe system crosses a saddle point instability, causing astress drop and particle rearrangements as the systemtransitions to a new energy minimum. The magnitude ofthe stress drop quantifies the size of the rearrangementevent.In AQS, the stress averaged over many such stressdrops gradually rises until about 6-7% strain, at whichpoint the systems yields. After the yielding point, theaverage stress remains constant as a function of strain.In this manuscript, we focus on the pre-yielding regime,where the infinite-dimensional mean-field equations aresolvable [8–12] [31], and where the response dependsstrongly on the initial preparation. Moreover, the lo-cal shear modulus µ , defined as the slope of the stress-strain curve along elastic branches, is significantly dif-ferent from the macroscopic coarse-grained shear modu-lus µ global , defined as the ratio of the average stress atyield to the average strain at yield. This observationis directly related to marginal stability [33], and can bequalitatively predicted from infinite dimensional analytictheory [8–12].To develop a more quantitative comparison betweenAQS and AQRD, as predicted in Eq. (7), we focus onthree metrics that quantify how AQS and AQRD samplephase space in the pre-yielding regime: (i) the distribu-tion of local shear/random-displacement moduli µ and ˜ µ along elastic branches, (ii) the distribution of (random)strain intervals ∆ γ and ∆˜ γ between stress drops, and (iii) the distribution of (random) stress drop magnitudes∆ σ and ∆˜ σ . We use ⟨ ∆˜ γ ⟩ and ⟨ ∆˜ σ ⟩ to denote quantitieswhich are explicitly averaged over all elastic branches inthe pre-yielding regime. FIG. 3.
System size and pressure dependence of land-scape statistics. a)
Local shear modulus µ , b) Strain dis-tance between rearrangements ∆˜ γ , and stress drops acrossrearrangements ∆˜ σ as a function of N p to show collapsewith system size and pressure. AQRD with completely un-correlated random fields is shown with closed circles, whileAQS data is shown with open circles. Error bars representthe middle 60% of the distribution and are only shown forAQRD for visual clarity, but are approximately the same forAQS. Colors represent system sizes N =
64 (red), 128, 256,512, 1024 (blue) in an even gradient.
Scaling of observables with system size and pressure—
Previous work has analyzed these statistics in AQS asa function of system size N and pressure p [27, 34, 35], as such data helps constrain continuum so-called ‘elasto-plastic’ models to predict features of avalanches in gran-ular matter. In addition, the size of a rearrangementprovides interesting information about the nonlinear fea-tures of the potential energy landscape, as it is one wayof quantifying how far the system has to travel from asaddle point to find a nearby local minimum. The size ofavalanches in AQS, quantified by the magnitude of thestress drops and other metrics, is known to exhibit power-law scaling with a large-scale cutoff, and the power lawhas different exponents on either side of the yielding tran-sition [34]. In the pre-yielding regime the average stressdrop is well defined, and changes in a systematic waywith system size and pressure. Previous work by some ofus [27] demonstrated that in AQS the average stress dropexhibits two regimes: a finite-size regime when N p ≪ N p ≫ ⟨ ∆ σ ⟩ ∼ pN , which is illustrated by the open symbols inFig. 3b.Therefore, we first study the statistics of stress dropsfor the simplest choice for the AQRD vector field ∣ c ⟩ –anuncorrelated random field (GCF with ξ = N and p , showing that precisely the samescaling is seen in AQRD. This highlights that the zero-pressure limit of the avalanche statistics under AQRD issingular, just as in AQS. Although the scaling is identicalthere is clearly a shift in the prefactors, which we returnto in the next section.In addition to the magnitude of the stress drops, thestrain between saddle points or rearrangements providesanother window into the statistical features of the com-plex potential energy landscape. Fig. 3c clearly showsthat the mean strain interval between rearrangementsscales as ⟨ ∆˜ γ ⟩ ∼ p / N in both AQS (open circles) andAQRD (closed circles). Additionally, we measure theaverage shear modulus between rearrangements, whichscales as ⟨ µ ⟩ ∼ p / for both AQS and AQRD as shownin Fig. 3a. Effects of spatially correlated forcing—
Although thescaling exponents of the previous section are precisely thesame under both AQS and AQRD dynamics, it is clearthat there is a systematic offset in the prefactors, despitethe fact that care was taken to ensure the definition ofeffective strain in each case is equivalent.To understand the origin of this difference, we vary thecorrelation length ξ of the normalized AQRD vector field ∣ c ⟩ measured in units of the smaller particle diameter anduse Gaussian correlated fields (GCF) and wave-like cor-related fields (WCF), as described in the methods sectionand illustrated in Fig. 2. For these analyses, system size N = φ = .
94 are fixed andknown to be far from the singular limit. These parame-
FIG. 4.
Collapse of landscape statistics with correlation length. a)
Probability distribution of the local effective moduli˜ µ and µ (inset) and the recentered z ˜ µ and z µ (rescaled by their standard deviations) in GCF systems with ξ = ξ = . µ of AQS (black). b) The average effective modulus decreases as a function of correlation lengthin both WCF and GCF ensembles. All curves approach the AQS value (black diamond), and dashed lines are best fits for˜ µ WCF with slope − . µ GCF with slope − . − − c) A comparison of F / (cid:96) computed directly via the variance of the field ∣ c ⟩ (black and gray lines)and the initial modulus ratio κ = ˜ µ / µ (magenta and green lines). d) Collapse of average stress-strain curves for GCF randomfields onto average AQS stress-strain curve using Eq. 9. Here ˜ σ avg and σ avg denote an average over configurations but not allelastic branches. Using κ ∗ = ⟨ ˜ µ ⟩/⟨ µ ⟩ and the data in b we are able to collapse the distributions of e) the effective strain interval∆˜ γ √ κ ∗ and f ) the effective avalanche size ∆˜ γ /√ κ ∗ , by appropriate scaling of the raw data (insets). Data shown is for GCF.The avalanche distribution agrees with the reported slope of − ters produce L x = L y = √ Nπ ( + . ) φ = .
3, where 1 and1 . ∣ c ⟩ for both GCF and WCF are shown inFig. 2. In each case, the random stress vs. random straincurves exhibit qualitatively similar features, with elasticbranches punctuated by stress drops. The overall magni-tude of the stress scale changes dramatically, where largerstresses are associated with smaller correlation lengths.In order to test the prediction of Eq. 7, we first in-vestigate the statistics of the local shear modulus, ˜ µ ,shown for the GCF data in the inset to Fig. 4a. TheGCF distributions shifted by the mean z ˜ µ and scaled bythe standard deviation do not collapse as shown in themain panel Fig. 4a. However, the average is well-definedfor both GCF and WCF data sets, and decreases withincreasing ξ (Fig. 4b). Specifically, both data sets areconsistent with ˜ µ being a power law function of ξ , andthe AQS data point shown by the black diamond falls onboth of the lines describing GCF and WCF data, respec-tively. This again confirms that shear is a special case ofa more generalized response to displacement fields.Next, we define a new variable, κ , as the initialrandom-displacement modulus ˜ µ normalized by the ini- tial shear modulus µ : κ ≡ ˜ µ / µ . We then explicitly testthe mean-field prediction for the shear modulus, Eq. 7: κ = ˜ µ / µ = F / (cid:96) . In order to compute these quantities inour simulation data, we follow the prescription of Eq. 6,taking F (cid:96) = dN c ∑ ⟨ i,j ⟩ ( c i − c j ) , (8)where N c is the total number of contacts, ⟨ i, j ⟩ denotescontacting neighbors, and we approximate (cid:96) as the aver-age distance between contacting particles. These quan-tities, calculated for both the Gaussian correlated fields F GCF and the wavelike correlated fields F WCF , are shownby the grey and black data points in Fig. 4c, respectively.We also plot the modulus ratio κ as a function of corre-lation length for both GCF(green) and WCF(magenta)simulations. Although this is a 2d system far from theinfinite-dimensional mean field case, the mean field pre-dictions are fairly close to the WCF data, and also cap-ture the general trend of the GCF data.However, the mean-field prediction is not in quantita-tive agreement so that F / (cid:96) ≠ κ , suggesting that in lowdimensions the dynamics cannot be reduced solely to the FIG. 5.
Average effective stress-strain curves in ultra-stable glasses . Stress-strain curves are shown with ξ = T init = .
062 (red), T init = . T init = . σ avg and σ avg denote anaverage over configurations but not all elastic branches. Thecurves are collapsed via Eq. 9. We see that for lower prepara-tion temperatures, there is a larger shear modulus and a morepronounced peak, in accordance with AQS simulations[36].Our predictions for the collapse agree well up to the yieldingpoint ( γ ≈ . variance of relative local strains. Nevertheless, a moregeneral prediction of the mean-field theory is that oncethe mechanical response at one value ξ is known, all oth-ers follow. Thus, one may expect that by using AQSas a reference state we can still collapse low-dimensionalsimulation data using κ : σ ∼ ˜ σ √ κ , γ ∼ ˜ γ √ κ. (9)For individual response curves, the proper value of κ isdefined at ˜ γ =
0, but to collapse landscape statistics,it suffices to use κ ∗ = ⟨ ˜ µ ⟩/⟨ µ ⟩ , which can be read fromFig. 4b, and which is explicitly averaged over all elas-tic branches in the pre-yielding regime, just as ⟨ ∆˜ σ ⟩ and ⟨ ∆˜ γ ⟩ . Physically, this means that the statistical featuresof the potential energy landscape are dominated by thescaling of the elastic moduli, i.e. by the curvature of thelandscape minima.Fig. 4d-f demonstrates that this mean-field predictionworks remarkably well: individual stress-strain curves,distributions of strain intervals between stress drops, andthe magnitude of stress drops all collapse when properlyscaled by κ as predicted by mean-field theory. In ad-dition, the collapsed avalanche data is clearly consistentwith the scaling of P ( ∆ σ ) ∼ / ∆ σ reported by Shang et.al [34]. This is another indication that bulk responses ofAQS and AQRD are controlled by the same physics. Effect of material preparation and stability—
Tothis point, we have investigated infinite-temperature-quenched jammed solids, which have a high degree of dis-order. Under AQS, such systems exhibit a ductile yield- ing transition where the pre-yielding regime transitionssmoothly to the post-yielding regime with no disconti-nuity in the stress. It is well-known that changing thematerial preparation protocol alters the disorder in theinitial configuration, and changes the yielding transition.Recent work using a new Swap Monte Carlo algorithmgenerates ultrastable glasses that are–on the contrary–extremely brittle, with large stress overshoots and dis-continuous stress drops at the yielding transition, anddata from such simulations strongly suggests that underAQS the yielding transition is in the Random Field IsingModel universality class [36, 37]. Although a full study ofthe nature of the yielding transition in AQRD is beyondthe scope of this work, we analyze the random-stress vs.random-strain curves using GCF under different prepa-ration protocols.The solid lines in Fig. 5a shows such curves for dif-ferent parent preparation temperatures, ranging from T init = . .
062 (brit-tle, ultrastable glass, high stability). The dashed curvescorrespond to the stress-strain response in AQS for thesame initial conditions. We observe that in AQRD, theglobal modulus increases as the stability increases, whichis similar to what is observed in AQS. In addition, there isclear stress overshoot (where the average stress increasesfar above its later steady state value) for the ultrastableglass, which is similar to what is seen for the yieldingtransition in AQS, although the yielding transition ismuch sharper in AQS. Taken together, these results high-light that the qualitative trends for how the yielding tran-sition depends on glass stability are similar in AQS andAQRD, and sets the stage for future work to study thestatistics and spatial structure of the yielding transitionin AQRD.
III. CONCLUSION
These results demonstrate that shear and randomforces perturb disordered solids in remarkably similarways. In particular, the nonlinear properties of the po-tential energy landscape traversed by AQS or AQRD dis-play identical scaling exponents. We discovered that theprefactors for these scaling laws, which generally charac-terize the stiffness of the material or the magnitude of thecurvature in the potential energy landscape, are a powerlaw function of the correlation length of the input fieldof displacements. The exponent ranges from − − IV. ACKNOWLEDGEMENTS
We thank Ludovic Berthier and Misaki Ozawa for dis-cussion and initial configurations of ultrastable glassesusing swap monte carlo. PM would like to thankMatthias Merkel for helpful discussions. EA wouldlike to thank Francesco Zamponi and Ada Altieri fordiscussions about the infinite-dimension mean-field re-sults. EA acknowledges support from the Swiss Na-tional Science Foundation by the SNSF Ambizione GrantPZ00P2 173962, from the European Research Council(ERC) under the European Union Horizon 2020 researchand innovation programme (grant agreement n. 723955- GlassUniversality), and from the Simons Founda-tion Grant ( [1] V. Narayan, S. Ramaswamy, and N. Menon, Science , 105 (2007).[2] M. E. Cates and J. Tailleur, Annu. Rev. Conden. Ma. P. , 219 (2015).[3] M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B.Liverpool, J. Prost, M. Rao, and R. A. Simha, Rev.Mod. Phys. , 1143 (2013).[4] S. Henkes, Y. Fily, and M. C. Marchetti, Phys. Rev. E , 040301 (2011).[5] S. Henkes, K. Kostanjevec, J. M. Collinson, R. Sknepnek,and E. Bertin, Nat. Commun. , 1405 (2020).[6] L. Berthier and J. Kurchan, Nat. Phys. , 310 (2013).[7] L. Berthier, E. Flenner, and G. Szamel, New J. Phys. , 125006 (2017).[8] C. Rainone, P. Urbani, H. Yoshino, and F. Zamponi,Phys. Rev. Lett. , 015701 (2015).[9] C. Rainone and P. Urbani, J. Stat. Mech. , 053302(2016).[10] P. Urbani and F. Zamponi, Phys. Rev. Lett. , 038001(2017). [11] A. Altieri and F. Zamponi, Phys. Rev. E , 032140(2019).[12] G. Parisi, P. Urbani, and F. Zamponi, Theory of SimpleGlasses: Exact Solutions in Infinite Dimensions (Cam-bridge University Press, 2020).[13] E. Agoritsas, T. Maimbourg, and F. Zamponi, J. Phys.A-Math. Theor. , 334001 (2019).[14] E. Agoritsas, T. Maimbourg, and F. Zamponi, J. Phys.A-Math. Theor. , 144002 (2019).[15] Q. Liao and N. Xu, Soft Matter , 853 (2018).[16] Y. Fily, S. Henkes, and M. C. Marchetti, Soft Matter , 2132 (2014).[17] P. Olsson and S. Teitel, Phys. Rev. Lett. , 178001(2007).[18] A. Ikeda, L. Berthier, and P. Sollich, Soft Matter , 7669(2013).[19] D. Bi, J. H. Lopez, J. M. Schwarz, and M. L. Manning,Nat. Phys. , 1074 (2015).[20] D. Bi, X. Yang, M. C. Marchetti, and M. L. Manning,Phys. Rev. X , 021011 (2016). [21] M. Merkel and M. L. Manning, New J. Phys. , 022002(2018).[22] N. Xu, V. Vitelli, A. J. Liu, and S. R. Nagel, Europhys.Lett. , 56001 (2010).[23] D. L. Malandro and D. J. Lacks, J. Chem. Phys. ,4593 (1999).[24] C. E. Maloney and A. Lematre, Phys. Rev. E , 016118(2006).[25] C. S. OHern, L. E. Silbert, A. J. Liu, and S. R. Nagel,Phys. Rev. E , 011306 (2003).[26] S. Dagois-Bohy, B. P. Tighe, J. Simon, S. Henkes, andM. van Hecke, Phys. Rev. Lett. , 095703 (2012).[27] P. Morse, S. Wijtmans, M. van Deen, M. van Hecke, andM. L. Manning, Phys. Rev. Res. , 023179 (2020).[28] E. Bitzek, P. Koskinen, F. Ghler, M. Moseler, andP. Gumbsch, Phys. Rev. Lett. , 170201 (2006).[29] P. K. Morse and E. I. Corwin, Soft Matter , 1248(2016).[30] P. Charbonneau, J. Kurchan, G. Parisi, P. Urbani, andF. Zamponi, Nat. Commun. , 3725 (2014).[31] E. Agoritsas (in preparation).[32] To compare the mean-field gap directly with the softspheres in our simulations, we can use the relationship h ij = − d(cid:15) ij ρ i + ρ j (cid:96) .[33] J. Lin and M. Wyart, Phys. Rev. X , 011005 (2016).[34] B. Shang, P. Guan, and J.-L. Barrat, Proc. Natl. Acad.Sci. U.S.A. , 86 (2020).[35] S. Franz and S. Spigler, Phys. Rev. E , 022139 (2017).[36] M. Ozawa, L. Berthier, G. Biroli, A. Rosso, and G. Tar-jus, Proc. Natl. Acad. Sci. U.S.A. , 6656 (2018).[37] M. Popovi, T. W. J. de Geus, and M. Wyart, Phys. Rev.E , 040901 (2018).[38] D. Richard, M. Ozawa, S. Patinet, E. Stanifer, B. Shang,S. A. Ridout, B. Xu, G. Zhang, P. K. Morse, J.-L. Barrat,L. Berthier, M. L. Falk, P. Guan, A. J. Liu, K. Martens,S. Sastry, D. Vandembroucq, E. Lerner, and M. L. Man-ning, arXiv:2003.11629 [cond-mat] (2020). Supplemental Material
I. GENERATING RANDOM VECTORS FORAQRDA. Wave-like correlated fields (WCF)
We generate the vector ∣ c ⟩ under the constraint thatit must be continuous across the boundary, and thus pe-riodic. In order to generalize AQS strain, we set the y -component of each c i to zero and let the x -componentdepend on the height of the particle c i = sin ( πy i ξ ) ˆ x (10)where ˆ x is the unit vector in the x -direction and ξ is thecorrelation length. In order for this to be periodic, wemust have ξ = L y / n with n ∈ Z . We note that AQS withLees-Edwards boundary conditions is simply ξ = L y witha phase shift (or simply a cosine). Once ∣ c ⟩ is determined,it is normalized so that ⟨ c ∣ c ⟩ = B. Gaussian correlated fields (GCF)
This section describes the process of generating a ran-dom Gaussian vector ∣ c ⟩ with finite spatial correlationlength ξ using Fourier transforms. This method is usedfor all GCF ensembles, except ξ = ∣ c ⟩ is drawn from a uniform distribution and then nor-malized to length 1.The two-dimensional system is a box of size L x × L y with periodic boundary conditions, allowing us to definewave vectors k nm = ( πnL x , πnL x ) for n, m ∈ Z . In practice,we truncate the Fourier sums at n, m = Q taking Q = ∣ c ⟩ , we need a correlated randomfield Ψ ( x ) which is Gaussian distributed with zeromean ⟨ Ψ ( x )⟩ = ⟨ Ψ ( x ) Ψ ( x ′ )⟩ = f (∣ x − x ′ ∣) . We enforce f ( x ) to be aGaussian function, whose explicit Fourier transform is˜ f (∣ k ∣) = exp [(− ∣ k ∣ ξ )] .First, we generate a set of uncorrelated random fields˜ ψ ( k ) with ⟨ ˜ ψ ( k )⟩ = ⟨ ˜ ψ ( k ) ψ ( k ′ )⟩ = π δ k , − k ′ ,where δ is the Kronecker function, and the factor 4 π comes from the Fourier transform convention. In prac-tice, for each wave vector of the truncated sum we gen-erate a random field ˜ ψ ( k nm ) = A ( k ) exp {( iB ( k ))} , with A ( k ) = A (− k ) normally distributed with zero mean andvariance π , and B ( k ) = − B (− k ) uniformly distributedon the interval [ , π ] .Secondly, we use the Fourier transform of thetarget correlator f ( x ) to construct a new field˜Ψ ( k ) = ˜ f (∣ k ∣) ˜ ψ ( k ) , whose Fourier transform is Ψ ( x ) = Q ∑ n,m = A nm e −∣ k nm ∣ ξ / cos ( B nm + k nm ⋅ x ) . (11)The random vectors c i are then defined as c i = Ψ ( x i ) using the initial positions x i . Once ∣ c ⟩ is determined, itis normalized so that ⟨ c ∣ c ⟩ = II. STRESS-CONTROLLED MEASUREMENTS
In addition to the strain-controlled definition of AQRDin the main text, we can also consider a stress-controlledversion. Instead of enforcing a displacement vector ∣ c ⟩ , weapply an external force ∣ F ext ⟩ = f ∣ c ⟩ on the system andwe measure the resulting strain, where again ⟨ c ∣ c ⟩ = F i = ∑ j ∈ ∂i F ij + f c i = , (12)where j ∈ ∂i indicates a sum over interparticle forcesacting on particle i . The sum of interparticle forces isparallel to ∣ c ⟩ , allowing closure with ∣ c ⟩ to give f = − ⟨ c ∣ F ⟩ .Following the arguments around Eq. (2) a normalizationfactor is necessary to compare the size of stresses in AQSleading to ˜ σ = − f L y √ N . (13)The random strain is then naturally defined as˜ γ = dUd ˜ σ = N ∑ i = ∂U∂ x i ⋅ d x i d ˜ σ . (14)where Eq. (12) ensures that ∂U∂ x i = ∑ j ∈ ∂i F ij = − f c i , and d x i d ˜ σ can be interpreted for small stresses as the displace-ment ∆ x i of particle i , subject to the applied stress ˜ σ ,leading to the natural relation d x i d ˜ σ = − ∆ x i f L y √ N . (15)Together these lead to the definition of random strainin the stress-controlled ensemble:˜ γ = ⟨ c ∣ ∆ x ⟩ L y √ N (16)Physically, this scalar product is simply the sum overall particles of their individual displacements along their2corresponding imposed direction c i . This is the natu-ral counterpart of the random stress definition given inEq. (4).In the linear response regime, Fig. S1 shows thatthe strain- and stress-controlled ensembles applied tothe same initial configuration give identical stress-straincurves, as expected. Beyond linear response, but still inthe pre-yielding regime (˜ γ < . FIG. S1. Comparison of strain controlled AQRD and stresscontrolled AQRD using the same initial configuration and thesame ∣ c ⟩ with N = φ = .
94 and a GCF correlationlength ξ ==