Probing relevant ingredients in mean-field approaches for the athermal rheology of yield stress materials

Although the notion of mechanical noise is expected to play a key role in the non-linear rheology of athermally sheared amorphous systems, its characterization has so far remained elusive. Here, we show using molecular dynamic simulations that in spite of the presence of strong spatio-temporal correlations in the system, the local stress exhibits normal diffusion under the effect of the mechanical noise in the finite driving regime. The diffusion constant appears to be proportional to the mean plastic activity. Our data suggests that the corresponding proportionality constant is density independent, and can be directly related to the specific form of the rheological flow curve, pointing the way to a generic way of modeling mechanical noise in mean-field equations.

Introduction -Amorphous materials, such as glasses or soft yield stress materials, exhibit a strongly heterogeneous response if driven by an external shear: fast particle rearrangements, the so-called shear transformations, take place in small regions while the rest of the material deforms elastically [1]. It is known that these plastic events induce long-range elastic deformations in the system [2,3]. In athermally sheared systems, where stress relaxation due thermal fluctuations can be neglected, like in the rheology of foams or gels, plastic deformation is controlled solely by the accumulation and interactions of these elementary events, that lead to avalanche type dynamics in the limit of slow driving [4,5].
Due to the long-range nature of the interaction between events, it is tempting to believe that the critical dynamics close to the yielding transition can be well described through mean-field considerations [6][7][8][9][10]. However, it has been shown in microscopic simulations [11] and using mesoscopic models [12,13] that athermally sheared systems can exhibit critical exponents that are different from the mean-field predictions, most probably due to the special form of the elastic kernel [13].
In this work we concentrate on the investigation of flow responses in a regime further away from the critical point, relevant for many rheological setups [14][15][16]. At higher shear rates and/or small enough system sizes we recover predictions made by mean-field models, in analogy with near mean-field critical points [17,18]. At the phenomenological level, there have been several propositions to model the mechanical noise [6,7]. We believe that the Hébraud-Lequeux (HL) model [8] and its further developments [9,19] provide so far the best selfconsistent mean-field description of mechanical noise in athermally sheared systems [20,21]. Within these models local stress diffusion is introduced with a diffusion coefficient proportional to the plastic activity, that can be * Francesco.Puosi@ens-lyon.fr measured experimentally [15]. In this letter we discuss the robustness of the HL picture on the basis of particle based simulations of a two dimensional athermally sheared amorphous solid. While our results support the modeling of mechanical noise via stress diffusion, other assumptions of the model, like the existence of a unique local critical stress and the complete stress release following a plastic event, are in clear contrast with microscopic observations. Our work is not intended as a definitive verdict on the validity of the HL model but rather as a starting point for further improvements of the model toward a more realistic description of the flow of disordered media.
Mean-field approach -The HL-model describes the state of a soft glassy material via the probability density P of local shear stresses σ in regions of mesoscopic size W while the material is sheared at rateγ. The time evolution of P is given by where θ(x) and δ(x) denote respectively the usual Heaviside and delta-distributions. The first term on the right hand side proportional to the stress gradient of the probability density ∂ σ P accounts for the linear elastic response.
The following term describes the loss in the probability density due to local yielding of overstressed regions above a critical stress σ c at a rate given by 1/τ . The corresponding gain term is given in the third expression on the right hand side, where the stress is set to zero after a yielding event with a rate given by the plastic activity rate Γ(t) = τ −1 dσP(σ, t)θ(|σ| − σ c ). The last term encompasses the stress changes created through other yielding events in a mean field manner, assuming that this mechanical noise can be approximated through a normal diffusion of the mesoscopic stresses. To describe this noise in an self-consistent manner, the HL approach proposes that its amplitude should be related to the rate of plastic activity through a coupling constantα  This model is known to exhibit a unique stationary state in the large time limit. Using appropriate units we can write the equations in dimensionless quantities, expressing stress related values in units of the local yield stress σ c , time quantities in units of τ , the shear rate in units of σ c /(G 0 τ ) and the stress diffusion coefficient in units of σ 2 c /τ , leaving only two independent adimensional model parameters that determine the flow behaviour, namely the dimensionless shear rate and the coupling constantα.
The rheological results in the small shear rate limit for this model are well studied [8,[21][22][23][24][25][26][27][28]. For a small enough coupling strengthα < 1/2 the HL model predicts a Herschel-Bulkley law of exponent 1/2, σ ≈ σ Y + Aγ 1/2 , for the dependence of the time averaged macroscopic stress in the steady state, given by σ = σP(σ)dσ. with the two constants σ Y (the dynamical yield stress) and A (the prefactor). These macroscopic constants of the model can only depend on the last free control parameter, namely the coupling constantα. This means that in dimensionless units for the prefactor and the dynamical yield stress we will obtain a universal curve A/ √ G 0 τ σ c versus σ Y /σ c parametrized throughα (see right panel of Fig. 1) [29].
In the following we aim at probing not only the above constant relation using microscopic dynamic simulations, but also to test the the underlying assumptions, most importantly the self-consistent description of the mechanical noise induced through the plastic activity. To test the above theory we need to measure all involved parameters that appear in the HL description. Parameters like the shear modulus G 0 are well defined quantities in the microscopic simulations and rather easy to measure. Further if stress fluctuations turn out to be diffusive, D HL has a well defined meaning. Other HL parameters like the local yield stress σ c , the phenomenological parameter τ and the rate of plastic activity Γ(t) are rather difficult to interpret within the microscopic picture. In the following we describe our attempt to relate the different parameters to measurable quantities in the microscopic dynamics.
Microscopic model -We have investigated a generic two-dimensional (2D) model of glass, consisting of a mixture of A and B particles, with N A = 10400 and N B = 5600, interacting via a Lennard-Jones potential V αβ (r) = 4ǫ αβ [(σ αβ /r) 12 − (σ αβ /r) 6 ] with α, β = A, B and r being the distance between two particles. The parameters σ αβ and ǫ αβ are chosen in order to prevent crystallization in 2D at low temperature [30]. The units of energy, length and mass are defined by ǫ AA , σ AA and m A ; the unit of time is given by The potential is truncated at r = r c = 2.5 for computational convenience and periodic boundary conditions are used. The equations of motion are integrated using the velocity Verlet algorithm with a time step δt = 0.005. The athermal limit is achieved by thermostating the system at zero temperature via a Langevin thermostat [31] with a damping coefficient ζ = 1.
By changing the volume of the system V, we explored states with different number density ρ = (N A + N B )/V . Glassy configurations were prepared by quenching equilibrated configurations at T = 1 to zero temperature with a fast cooling rate. Simple shear is set at a ratė γ by deforming the box dimensions and remapping the particle positions. In order to characterize the plastic activity of the system we consider the D 2 min quantity [7]. For a given particle i, D 2 min is defined as the minimum over all possible linear deformation tensors ǫ of: where the index j runs over all the neighbors of the reference particle i and I is the identity matrix. To temporally resolve the plastic events, we set the time lag to the value δt = 4.
Macroscopic flow curve -In Figure 1(a) we show the dependence of the macroscopic shear stress σ xy on the applied shear rateγ. The flow curves are well described by the Herschel-Bulkley (HB) law, σ xy = σ Y + Aγ n , with an exponent n ≈ 0.55 (not shown here) which seems not to depend on the density. Fixing the exponent n to the value 0.5, the one predicted by the HL model in the case ofα < 1/2, gives indistinguishable results (see Fig. 1(a)). Although other works on sheared disordered material in two and three dimensions report similar values for the flow curve exponent [32,33], recent works in the literature seem to suggest that a proper finite size scaling analysis close to the yielding transition reveals different critical dynamical exponents [12,13,34]. We would like to insist here on the fact, that our study does not aim at measuring the critical exponents of the transition; instead we rather test the consistency of the assumptions made in the HL approach in a parameter regime that fits well the model predictions.
Size of an elementary plastic region -Here we describe the procedure we followed to convert microscopic observations into mesoscopic ones. First, we define a microscopic criterion for active particles, i.e., particles performing plastic rearrangements, based on the D 2 min quantity: a particle i is active at a given time t if the corresponding D 2 min (i, t) is larger than a threshold value that we fix equal to 0.1. In Figure 2(a) we compare a typical stress-strain curve and the corresponding evolution of the total number N pl of active particles. The correlation between the stress drops and the peaks in N pl suggests that the actual definition, in particular the threshold we introduced, is reasonable. Then a high-resolution discretization of the system is performed by dividing it into k × k cubic blocks of size w 0 = 2. A small block is considered as active if it contains at least one active particle.
The mean size of a plastic events l pl can be estimated by a cluster analysis of the active blocks in the configurations explored by the system, namely l pl = A pl

1/2
where A pl is the volume of a plastic cluster. We obtain l pl ≈ 6, with no relevant dependence on the density. This value is in accord with previous works reporting plastic regions with a size of a few particle diameters [32,35]. Moreover, in the rangeγ 10 −4 , l pl does not exhibit a strong dependence on the shear rate. Next we implement the coarse-graining of microscopic simulations on the scale of individual plastic events. The simulation box is divided into M × M cubic blocks with M chosen in order to have W ≈ l pl . The local stress σ m αβ of a cube m is defined as: whereρ m is the number density in the cube m and the summation of i is performed over the particles in the cube. The macroscopic stress tensor σ αβ is obtained by the summation of σ m αβ over all the blocks. Stress diffusion -We define the coarse mean-square stress difference as: where σ m xy (t) is the stress in a given block at time t and the last term in Eq. 4 account for the stress increase due to the elastic deformation of the system, being G the macroscopic shear modulus. In Figure 2, the short time behavior of (∆σ m xy ) 2 (t) is shown for a finite shear rateγ = 10 −6 , approaching the quasistatic limit. We observe that at short times, (∆σ m xy ) 2 (t) increases linearly with time. We define the stress diffusion coefficient as D σ = (∆σ m xy ) 2 (t) /t. In the inset of Fig. 2 we show the dependence of D σ on the density of the system.
Duration of a plastic event -In Fig. 2(c) we show the D 2 min autocorrelation function for active particles C p = D 2 min (t 0 )D 2 min (t 0 + t) / (D 2 min (t 0 )) 2 as a function of time. We observe that C p decays exponentially with a characteristic time τ p , that depends weakly on density (see the inset of Fig. 2) and that is close to the damping time τ d = Γ −1 of the Langevin thermostat. Local yield stress -In this section, we present a method which allows us to calculate a local critical stress, i.e., the stress limit before a plastic rearrangement occurs locally. For this purpose, we adopt the "frozen matrix" method [36][37][38]. The system is frozen except for a target region, i.e., a mesoscopic block. When a deformation is applied, the frozen region can only deform affinely whereas the target block is allowed to relax non-affinely. Hence, the local strain ǫ m of the target block coincides with applied global strain while the local stress σ m αβ is given by Eq. 3. The system is submitted to a simple shear deformation following a quasi-static protocol. The strain is increased with a step δǫ; after each increment the system is relaxed by the steepest descendent method with the restriction that only the particles in the target region can move. Here we used δǫ = 10 −5 .
For small deformations the target region behaves elastically and the stress increases linearly, with a slope controlled by the local shear modulus [38]. As the strain increases, the elastic behavior goes on until the local yield stress is reached and a plastic rearrangement takes places. This is indicated by a stress drop in the stressstrain curve. For a given block, we define as the local yield stress σ m c the value of the local stress at the first maximum.
In Fig. 2(d) we show the distribution of σ m c for the different values of the density. First, the fact the local yield stress is distributed is in clear contrast with the HL assumption of a unique critical stress σ c . The exposed dis- With the frozen matrix method we estimated also the stress release following a local yielding event. We observe that the plastic event only partially relaxes the accumulated stress (not shown) in contrast with the assumption of the HL model of a complete relaxation. The fraction of the relaxed stress seems not to depend on density being σ m c / ∆σ m ≈ 5. In Ref. [38] it was shown that the mean value of the shear modulus obtained with the frozen matrix method depends on the target region size W and that it converges, from higher values, to the macroscopic modulus as W increases. This is due to the frozen environment which reduces the non-affine motion of the particles in the target region. For W ≈ 6 the discrepancy is less than 50%, a significant error. We are aware that the estimate of σ m c may be affected by this effect. Also we have to be careful in the interpretation of the result, because by using the frozen matrix method we measure the yield stress distribution, obtained from strained configurations. It would be interesting to try to infer the inherent yield stress distribution from our measurements [21], but this is beyond the scope of this study and left for further investigations.
Robustness of HL model -To probe the consistency of the HL predictions we examine first the theoretical constant correlation between the dimensionless yield stress σ Y /σ c and prefactor A/(G 0 σ c τ ) 1/2 in the flow curve prediction of the HL model. Yield stress and prefactor, as well as the shear modulus G 0 can be easily accessed in simulations, the formers from a HB fit of the flow curve, the latter as the slope of the stress-strain curve in the very first elastic regime. We approximate the critical stress σ c as the mean value σ m c of the local yield stress obtained with the frozen matrix method. As discussed in a recent work on the HL model, the phenomenological parameter τ can be interpreted as the duration of a plastic event in the low shear rate limit [21]. Hence, we will approximate τ by the value of τ p obtained through the measurement  of the D 2 min two-time autocorrelation function. In Fig. 1(b) we compare the dimensionless yield stress and prefactor from simulations with the theoretical correlation curve. First, we observe that all the data belonging to different density values collapse onto a single point very close to the theoretical curve. This means that although the flow curves are sensitive to the density, we find that the coupling constantα flow appears density independent, which might point to a universal coupling constant for this system that probably only depends on the specific form the the elastic propagator [3], in agreement with the KEP model predictions [9].
Although this first comparison appears very convincing we will now probe the second HL prediction, namely the self-consistent description of the noise. The value ofα flow obtained from the flow curve for the different densities considered are given in table I and yield a value of approximately 1/3. This is to be compared to the coupling constant obtained from its original definition as the proportionality factor between the diffusion coefficient and the plastic activity,α micro = D σ /(σ 2 c Γ) (see Eq. (2)), that we obtain from simulations with Γ from the cluster analysis. It turns out thatα micro is one magnitude smaller than the value calculated from the flow curve, indicating an inconsistency of the HL model. However, we note that both measurements are in good approximation density independent. This points toward the fact that the difference in the measurement ofα most probably results from oversimplifications in the HL dynamics.
It has been shown that the introduction of a local yield stress distribution does not alter the predicted HL exponent in the Herschel Bulkley approximation [21]. And we expect only slight changes in the values of the dimensionless yield stress and prefactor that cannot be responsible for the large difference in the measurement of the coupling constant. We expect that the partial relaxation of the stress after a yielding event, σ m c / ∆σ m ≈ 5, has a much stronger corrective effect. If we assume the diffusion coefficient to be proportional to the square of a typical local stress jump devided by a typical time scale given by the inverse plasticity rate, we obtain for the stress diffusion an altered expression Thus we obtain a new estimation for the coupling constantα ′ micro = ( σ m c / ∆σ m ) 2 α micro , which compares now much better with the coupling constant measured using the macroscopic flow curve (see table I). This suggests that the simplification of setting the stress to zero after a yielding event in the HL model has important consequences on the tested relation. Revisiting this rule by changing the gain term in the evolution equation (1) together with a more realistic yield stress distribution seems to us a very promising route to reach a more realistic mean field modeling of athermally sheared amorphous systems.
The most intriguing result of this study remains however the finding, that the coupling constant appears den-sity independent and as such might be a good candidate to classify the yielding properties of this system. It would be very interesting to be able to link the value of the coupling constant to mesoscopic measurements of the propagator as proposed in the KEP model and to perform similar analysis for systems with other dynamics to probe the generality of the concept.