Probing relevant ingredients in mean-field approaches

for the athermal rheology of yield stress materials

Francesco Puosi,∗a,b,c Julien Olivier,d and Kirsten Martensb,c
arXiv:1501.04574v3 [cond-mat.soft] 26 Aug 2015

1 Introduction
The theoretical understanding of the yielding transition in ather-
mally driven disordered systems is a highly challenging problem
and no consensus has been established even on the basic ingredi-
ents that should underlie coarse grained descriptions of the non-
linear rheological response of yield stress materials 1–6 .
The only commonly accepted point of view is that disordered
materials, such as glasses or soft matter, exhibit a strongly hetero-
geneous dynamics when 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 elasti- Fig. 1 Schematic view of the coarse grained picture: The particle
simulation box is devided into smaller mesoscopic parts that can hold
cally 7 . These plastic events induce long-range elastic deforma- exactly one plastic event. The plastic activity in the surrounding of a
tions in the system leading to complex correlations of the yielding mesoscopic region (see graphical interpretation on the right) leads to
regions in form of plastic avalanches 8–14 . stress fluctuations around the mean value of the stress (dashed line in
One of the major concerns remains to know whether despite the graph on the right), which is controlled by the external forcing.
this dynamical complexity it is sensible to describe the yielding
dynamics within mean-field descriptions 1–3,15,16 . We tackle this
key question using particle based simulations, concentrating on In this work we show that the dynamically created noise in
flow responses in a regime that is relevant for many rheological sheared amorphous systems can be encompassed through a nor-
setups 17? –26 . We find that at high enough shear rates and/or mal diffusion in the local stresses 33 (for a schematic view see
small enough system sizes we recover a dynamics well described Fig. 1). In agreement with former works 5,6 , we consider that
by mean-field considerations 27 , similar to near mean-field critical local stress diffusion is acting as noise in the tilt of the local po-
points in equilibrium phase transitions 28,29 . tential energy landscape with an amplitude proportional to the
An important feature that found its way into several theo- plastic activity, an observable that can be measured experimen-
ries, e.g. the Soft Glassy Rheology (SGR) and the Shear Trans- tally 21,24 .
formation Zone Theory (STZ), is to describe yielding in a gen- We analyse within molecular dynamics simulations the cou-
eralized thermodynamic description with an effective tempera- pling strength between the mechanical noise and the plastic ac-
ture 1,3 . But the analogy with thermally activated events proposed tivity. Within our simulations we successfully relate this cou-
by the SGR theory has recently been questioned for athermal rhe- pling strength, a dimensionless and density independent quan-
ology 5,6 . Also it has been shown that the thermodynamic inter- tity, in a consistent manner to the flow response without any
pretation of the effective temperature is problematic in the ather- further parameters. This finding leads us to reconsider a spec-
mal regime 30,31 and that the STZ theory is not able to predict the ulative mean-field scenario put forward a long time ago in the
non-linear response 32 , expected in the small driving limit 27 . so-called Hébraud-Lequeux (HL) model 2 , that was at the basis
of important further developments 15,22,34 , providing so far one
of the best self-consistent mean-field description of mechanical
∗
Corresponding author: francesco.puosi@ens-lyon.fr noise in athermally sheared disordered systems.
a
Laboratoire de Physique de l’École Normale Supérieure de Lyon, Université de Lyon,
CNRS, 46 Allée d’Italie, 69364 Lyon cédex 07, France. 2 Mean-field approach
b
Université Grenoble Alpes, LIPHY, F-38000 Grenoble, France.
c Within the framework of the HL-model we can establish a link be-
CNRS, LIPHY, F-38000 Grenoble, France.
d
Aix Marseille Université, Centre Mathématiques et Informatique, Technôpole Château tween the dynamical yield stress and the prefactor of a Herschel
Gombert, 39 Rue F Joliot Curie, F-13453 Marseille 13, France. Bulkley type power law fit of the rheological curves of athermally

1–8 | 1
sheared yield stress materials. This relation should, according
to this theory, solely be determined by the coupling strength be- 0.3 density
tween mechanical noise and the rate of activity. In the following 1.149

σxy
1 1.189
we review briefly the main assumptions and results of the mean- 1.253
1.322
field approach 2 . The model describes the state of a soft glassy ma-

σxy/σc
terial via the probability density P of local shear stresses σ = σxy -6 -4 -2
10 10 10
in regions of mesoscopic size W while the material is sheared at .
γ
rate γ̇. The time evolution of P is given by

∂t P(σ ,t) = −G0 γ̇(t)∂σ P − τ1 θ (|σ | − σc )P 0.2

+Γ(t)δ (σ ) + DHL (t)∂σ2 P (1) -4 -2

10 10
.
where θ (x) and δ (x) denote respectively the usual Heaviside and γ (G0τ/σc)
delta-distributions. The first term on the right hand side pro-
portional to the stress gradient of the probability density ∂σ P Fig. 2 Flow curves: Rescaled macroscopic shear stress σxy /σc as a
accounts for the linear elastic response. The following term de- function of the rescaled applied shear rate γ̇(G0 τ0 /σc ) for different
scribes the loss in the probability density due to local yielding of densities (see text). The dashed line is a fit with the Herschel-Bulkley
form σxy = σY + Aγ̇ 1/2 . The inset shows the same data before rescaling,
overstressed regions above a critical stress σc at a rate given by
dashed lines display a guide to the eye with the same Herschel-Bulkley
1/τ. It has been argued that the phenomenological parameter τ fitting form.
can be interpreted as the duration of a plastic event in the low
shear rate limit 6 . The corresponding gain term is given in the
third expression on the right hand side, where the stress is set to less shear rate and coupling constant α̃.
zero after a yielding event with a rate given by the plastic activity The rheological results in the small shear rate limit for this
rate model are well studied 2,6,35–41 . For small enough coupling
1
Z
Γ(t) = P(σ ,t)dσ . (2) strength α̃ < 1/2 the HL model predicts a Herschel-Bulkley flow
τ |σ |>σc
behaviour of exponent 1/2, hσ i ≈ σY + Aγ̇ 1/2 , for the time aver-
aged macroscopic stress in the steady state hσ i = σ P(σ )dσ ,
R
The last term encompasses the stress changes created through
other yielding events in a mean-field manner, assuming that this with the two constants σY (the dynamical yield stress) and A (the
mechanical noise can be approximated through a normal diffu- prefactor). These macroscopic constants of the model can only
sion of the mesoscopic stresses. To describe this noise in an self- depend on the last free control parameter, namely the coupling
consistent manner, the HL approach proposes that its amplitude constant α̃. This means that we will obtain a universal curve∗
√
should be related to the rate of plastic activity through a dimen- for the rescaled quantities A/ G0 τσc versus σY /σc parametrized
sionless coupling constant α̃ through α̃ (see Fig. 3).
The physical interpretation of this coupling constant could be a
DHL (t) = α̃σc2 Γ(t) . (3) mechanical fragility 2 that would depend on the details of the mi-
croscopic interactions and thus should be material dependent. In
This last relation introduces a non-linearity into Eq. (1), since a later work by Bocquet et al., the authors use a spatial approach
the rate of plastic activity itself depends on the density probability to the problem that allows for the derivation of an expression
of the mesoscopic stresses. It is this coupling that renders the for the coupling constant α̃ if one assumes a decorrelation of the
problem non-trivial and yields interesting results regarding the plastic event dynamics. Within this approximation it is possible
behaviour of macroscopic quantities. to express α̃ as a function of the specific form of the elastic stress
This model is known to exhibit a unique stationary state for propagator 15 . This suggests that there should be classes of ma-
a finite shear rate in the large time limit, where the probability terials with similar elastic responses that not only share the same
density for the mesoscopic stresses becomes time independent. non-linear exponent in the flow curve, but also comparable rela-
To determine the time averaged macroscopic stress in the steady tions between the dynamical yield stress and the prefactor in the
state one averages over the mesoscopic stresses weighted by the flow curve.
corresponding steady state probability density In the following we aim at probing not only the above con-
Z stant relation using microscopic dynamic simulations, but also to
hσ i = σ P(σ )dσ . (4) test the the underlying assumptions, most importantly the self-
consistent description of the mechanical noise induced through
Using appropriate units we can write the equations in dimension- the plastic activity. To test the above theory we need to measure
less quantities, expressing stress related values in units of the lo- all involved parameters that appear in the HL description.
cal yield stress σc , time quantities in units of τ, the shear rate in Quantities and parameters like the local stress and the shear
units of σc /(G0 τ) and the stress diffusion coefficient in units of
σc2 /τ, leaving only two independent dimensionless model param-
eters that determine the flow behaviour, namely the dimension-
∗ For the analytical derivation of this constant relation see the Appendix A

2| 1–8
 4 ing the volume V of the system. Glassy configurations were pre-
(a)
A (b)
4
pared by quenching equilibrated configurations at T = 1 to zero
temperature with a fast cooling rate. Simple shear is set at a rate
3 density 3
γ̇ by deforming the box dimensions and remapping the particle
1/2

1.149
2
A/(G0 τ σc)

1.189 positions. The quenching protocol has virtually no effect since

1.253 0 1 σ 2
y
2 1.322 we focus on the steady state shear (total imposed deformation
0.4 (c)

1/2
∆γ > 20%).

A/(G0 τ σc)
0.3
0.2
In order to characterize the plastic activity of the system we
1
0.1 consider the D2min quantity 3 . For a given particle i, D2min is defined
0.2 as the minimum over all possible linear deformation tensors ε of:
σy/σc
2
0 D2 (i,t, δt) = ∑ ri j (t + δt) − (I + ε) · ri j (t)

0 0.1 0.2 0.3 0.4 0.5 (5)
j
σy/σc
where the index j runs over all the neighbors of the reference
particle i and I is the identity matrix. We set the time lag to
Fig. 3 HL predictions vs. MD simulations: (a) The points are the rescaled
Herschel Bulkley prefactor A as a function of the rescaled dynamical δt = 4. This value is a compromise between having a good sig-
yield stress σY for four different densities, error bars are estimated from nal, i.e. large irreversible displacements, and being able to resolve
the measurements of the different quantities. The full line corresponds individual plastic events.
to the HL prediction. (b) This inset shows the data from the simulation
fits before the rescaling. (c) This inset shows a zoom into the main
panel, to display better the proximity to the theoretical curve obtained 3.1 Macroscopic flow curve
from a very simplified picture. In Fig. 2 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 expo-
modulus G0 are well defined quantities in the microscopic sim-
nent n ≈ 0.55 (not shown here) which seems not to depend on the
ulations and rather easy to measure. Also, if stress fluctuations
density. Fixing the exponent n to the value 0.5, the one predicted
turn out to be diffusive, DHL has a well defined meaning. How-
by the HL model in the case of α̃ < 1/2, gives indistinguishable
ever, other HL parameters like the local yield stress σc , the phe-
results. Although other works on sheared disordered material in
nomenological parameter τ and the rate of plastic activity Γ(t) are
two and three dimensions report similar values for the flow curve
rather difficult to interpret within the microscopic picture. In the
exponent 45,46 , recent works in the literature seem to suggest that
following we describe our attempt to relate the different parame-
a proper finite size scaling analysis close to the yielding transi-
ters to measurable quantities in the microscopic dynamics.
tion reveals different critical dynamical exponents 13,14,27,44 . We
3 Microscopic model would like to insist here on the fact, that our study does not aim
at measuring the critical exponents of the transition; instead we
We have investigated a generic two-dimensional (2D) model of
rather test the consistency of the assumptions made in the HL ap-
a glass, consisting of a mixture of A and B particles, with NA =
proach in a parameter regime that fits well the model predictions.
10400 and NB = 5600, interacting via a Lennard-Jones potential
Vαβ (r) = 4εαβ [(σαβ /r)12 − (σαβ /r)6 ] with α, β = A, B and r be-
ing the distance between two particles. The parameters εAA , σAA 3.2 Size of an elementary plastic region
and mA define the units of energy, length and mass; the unit of Here we describe the procedure we followed to convert the micro-
p
time is given by τ0 = σAA (mA /εAA ). We set εAA = 1.0, εAB = 1.5, scopic simulations into a mesoscopic description. First, we denote
εBB = 0.5, σAA = 1.0, σAB = 0.8 and σBB = 0.88 and mA = mB = 1. a particle i as active, i.e. performing a plastic rearrangement, at a
This choice is known to prevent crystallization in 2D at low tem- given time t if the corresponding D2min (i,t) is larger than a thresh-
perature 42 . The potential is truncated at r = rc = 2.5 for compu- old value that we fix equal to 0.1. In Fig. 4(a) we compare a typi-
tational convenience and periodic boundary conditions are used. cal stress-strain curve and the corresponding evolution of the total
The equations of motion are integrated using the velocity Ver- number N pl of active particles. The correlation between the stress
let algorithm with a time step δt = 0.005. The athermal limit is drops and the peaks in N pl suggests that the actual definition is
achieved by thermostating the system at zero temperature via a reasonable. Then a high-resolution discretization of the system is
Langevin thermostat 43 with a damping coefficient ζ = 1 which performed by dividing it into k × k square blocks of length w0 = 2.
corresponds to a strongly overdamped condition for the dynam- A small block is considered as active if it contains at least one
ics 44 . This model, which has been widely studied in different active particle.
versions 3,10,44,45 , is usually considered appropriate for colloidal The mean size, i.e. the mean linear extension, of a plastic
and other soft glasses. This choice is motived by the purpose of event hl pl i can be estimated by a cluster analysis of the spatial
investigating the general aspects of the rheology of athermal yield arrangement of the active blocks in the configurations explored
stress materials. by the system. We employed a modified version of the Hoshen-
To investigate different athermal flow responses, we explore Kopelman algorithm ? in order to account for periodic boundary
states with different number density ρ = (NA + NB )/V by chang- conditions. If we assume hl pl i = hA pl i1/2 where A pl is the area of

1–8 | 3
a plastic cluster, we obtain hl pl i ≈ 6, with no relevant dependence test this idea we define the coarse mean-square stress difference
on the density and shear rate (in the range γ̇ 6 10−4 ). The value as:
m 2
is in accord with previous works reporting plastic regions with h(∆σxy m
) (t)i = h(σxy m
(t0 + t) − σxy (t0 ) − G0 γ̇t)2 i (7)
a size of a few particle diameters 45,47 . Furthermore, this agree- m (t) is the stress in a given block at time t and the last
where σxy
ment justifies the criterion we employed to define active particles.
term in Eq. (7) accounts for the stress increase due to the elas-
Indeed, if we improve the resolution in the analysis of the plas-
tic deformation of the system, being G0 the macroscopic shear
tic activity by decreasing the threshold on D2min by a factor 10, m )2 (t)i
modulus. In Fig. 4(b), the short time behavior of h(∆σxy
while the number of active particles increase by a factor 4, the −6
is shown for a finite shear rate γ̇ = 10 , approaching the qua-
mean extension of a plastic event hl pl i is reduced by half, suggest- m )2 (t)i in-
sistatic limit. We observe that at short times, h(∆σxy
ing that single particle rearrangements are erroneously taken into
creases linearly with time. We define the stress diffusion coef-
account. m )2 (t)i/t. In the inset of Fig. 4(b) we show
ficient as Dσ = h(∆σxy
Next we implement the coarse-graining of microscopic simula-
the dependence of Dσ on the density of the system.
tions on the scale of individual plastic events. The simulation box
To define a duration of a plastic event we analyse the two-
is divided into M × M square blocks with M chosen in order to
m of a block m is defined
time autocorrelation function of the D2min quantity for active
have W ≈ hl pl i. The local shear stress σxy
particles. In Fig. 5(a) we show this correlation function C p =
as:
N ij ij hD2min (t0 )D2min (t0 + t)i/h(D2min (t0 ))2 i as a function of time. We ob-
m 1 1 ∂V (ri j ) rx ry
σxy = − 2 ∑ mvi,x vi,y + ∑ ∑ (6) serve that C p decays exponentially with a characteristic time τ p ,
W i∈m 2W 2 i∈m j=1 ∂ ri j ri j
that depends weakly on density (see the inset of Fig. 5(a)) and
where vi,x and vi,y are the x and y components of the velocity of that is close to the damping time τd = ξ −1 of the Langevin ther-
the particle i, ri j is the distance between the particles i and j and mostat. We choose to interpret this decorrelation time as the typ-
the summation of i is performed over the particles in the block. ical duration of a plastic event entering the HL model description.
The macroscopic stress tensor σxy is obtained by the summation
m over all the blocks.
of σxy 3.4 Local yield stress
In this section, we present a method which allows us to calculate
(a) a local critical stress, i.e., the stress limit before a plastic rear-
300
G0 rangement occurs locally. For this purpose, we adopt the “frozen
200 matrix" method 48–50 . The system is frozen except for a target re-
σxy

0.5 Npl
gion, i.e., a mesoscopic block, and it’s subjected to a simple shear
100 deformation, with a quasi-static protocol. The frozen region can
only deform affinely whereas the target block is allowed to relax
0.45 non-affinely. For small deformations the target region behaves
0 0.002 0.004
∆γ elastically and the stress increases linearly, with a slope controlled
10
0
by the local shear modulus 50 . As the strain increases, the elastic
(b)
~t behavior goes on until the local yield stress is reached and a plas-
<(∆σxy) (t)>

-2
10 tic rearrangement takes places. This is indicated by a stress drop
m 2

Dσ
-4 in the stress-strain curve. For a given block, we define as the local
10 density
-5
8×10
1.149
yield stress σcm the value of the local stress at the first maximum.
-5
-6 4×10
10 1.189
1.253
In Fig. 5(b) we show the distribution of σcm for the different
0
-8
1.322 1.2 ρ 1.3 values of the density. First, the fact the local yield stress is dis-
10
10
1
10
2 3
10 10
4 tributed is in clear contrast with the HL assumption of a unique
t critical stress σc . The exposed distributions are well described by
a simple Gaussian forms. The mean values hσcm i and the variance
Fig. 4 Activity and stress diffusion: (a) The full line represent a part of hδ σcm i of the distributions are shown in the inset of the figure. As
the macroscopic stress-strain curve in the steady state regime for a the density increases, hσcm i increases, due to the enhancement of
density ρ = 1.149 and a shear rate γ̇ = 10−6 . The average slope of the the repulsive interactions between particles.
elastic parts on the curve yields G0 . Open symbols: corresponding
evolution of the number of active particles N pl in the system. (b) Coarse With the frozen matrix method we estimated also the stress
mean-square stress difference h(∆σxy m )2 (t)i as a function of time for release following a local yielding event. We observe that the plas-
different values of the density. Note the linear behaviour at short times. tic event only partially relaxes the accumulated stress in contrast
Inset: stress diffusion coefficient Dσ = h(∆σxy m )2 (t)i/t as a function of
with the assumption of the HL model of a complete relaxation.
density.
The fraction of the relaxed stress seems not to depend on density
being h∆σ m i/hσcm i ≈ 0.2 (see inset of Fig. 5).
In a former work 50 it was shown that the mean value of the
3.3 Stress diffusion and duration of a plastic event shear modulus obtained with the frozen matrix method depends
The above introduced mean-field model assumes local stress fluc- on the target region size W and that it converges, from higher
tuations, obeying a normal diffusion process in stress space. To values, to the macroscopic modulus as W increases. This is due to

4| 1–8
 1 tion of the mechanical noise, thus predict a diffusion coefficient
(a)
1.3 proportional to the square of this typical local stress jump ∆σHL
densities
divided by a typical time scale, given by the inverse plasticity rate
Cp(t)

1.149 1.2
1.189 τp Γ(t).
0.5
1.253
1.322 (a) Our data analysis reveals that a homogeneous yield stress is
1.1
of course not verified in a disordered system, where one expects
1.2 1.3
ρ a distribution of yield stresses (see Fig. 5b). However, it has been
0 recently shown that the existence of a local yield stress distribu-
0 5 10 15
t tion in the HL dynamics does not strongly alter the predictions for
the flow behaviour 6 .
(b) 10 2
0.6 (b) We tested as well the second assumption of a normal diffu-

<∆σ >/<σc >

m
<δσc >
<∆σ >
<σc >

sion of the mesoscopic stresses, and we find that within the range

m
m
P(σc )

m
m

of shear rates that we consider, the measurement of the average

0 0
0.3 1.2 ρ 1.4 1.1 1.2 ρ 1.3 mean-square stress differences indeed allows for the determina-
tion of a well-defined diffusion coefficient. The HL coupling con-
HL = D /(σ 2 Γ) (see Eq. (3)) with D from the results
stant α̃micro σ c σ
0 in Fig. 4(b) and Γ from the cluster analysis turns out to be very
0 5 10 15 HL of order 10−2 for all considered densities.
m small, with α̃micro
σc (c) However, our study suggests that the partial relaxation of
the local stresses after yielding, h∆σ m i/hσcm i ≈ 0.2 (instead of the
Fig. 5 Duration of events and local yield stress: (a) D2min autocorrelation assumed total relaxation) introduces an important corrective fac-
function for active particles C p = hD2min (t0 )D2min (t0 + t)i/h(D2min (t0 ))2 i as a tor for the coupling constant. Since the typical stress jump in the
function of time. Lines are fitting curves exp(−t/τ p ). Inset: dependence diffusion process is now reduced, we obtain for the stress diffu-
of τ p on the density. (b) Distribution of the local yield stress p(σYm ) for
sion coefficient an altered expression
different values of the density. Full lines are Gaussian distributions fits.
Left inset: Mean value hσcm i and variance hδ σcm i as a function of density.
Dashed lines are power-law guides. Right inset: Absolute local stress Dσ = α̃micro (∆σ m )2 Γ. (8)
drop h∆σ m i and the relative local stress drop h∆σ m i/hσcm i as a function
of density. Dashed lines are power-law and constant value guides to the We report the values obtained for this new definition of the cou-
eye. Symbol code is the same as in panel (a). pling constant in Table 1. It turns out that these values are in good
approximation density independent and yield an average value of
about α̃ ≈ 0.26. Thus the coupling constant displays a reasonable
the frozen environment which reduces the non-affine motion of value smaller than 0.5 for which the HL model predicts 6 a well-
the particles in the target region. For W ≈ 6 the discrepancy is up established Herschel Bulkley regime in the flow curve with an
to 50%, a significant error. We are aware that the estimate of hσcm i exponent of 1/2 .
may be affected by this effect. Also we have to be careful in the And indeed our flow curves from the microscopic simulations
interpretation of the result, because by using the frozen matrix can be well fitted with an Herschel Bulkley expression σxy =
method we measure the yield stress distribution, obtained from σy + Aγ̇ 1/2 (see inset of Fig. 2). In Fig. 2 we plot our data in
strained configurations. It would be interesting to try to infer the dimensionless units, rescaled as suggested by the HL model, us-
inherent yield stress distribution from our measurements 6 , but ing dimensionless quantities σxy /σc and γ̇(G0 τ/σc ). The shear
this is beyond the scope of this study and left for further investi- modulus G0 can be easily accessed in simulations as the slope
gations. of the stress-strain curve in the elastic regime. We approximate
the critical stress σc as the mean value hσcm i of the local yield
4 Robustness of the HL model stress obtained with the frozen matrix method and we associate
In the following we aim at probing the robustness of the HL model the phenomenological parameter τ with the duration of a plastic
by comparing our data from the microscopic simulations to the event in the low shear rate limit 6 . Hence, we will approximate τ
various assumptions and predictions of the mean-field descrip- by the value of τ p obtained through the measurement of the D2min
tion. two-time autocorrelation function. We observe a collapse of the
Let us first recall some of the basic assumptions made to write flow curves for different densities. This collapse is not perfect,
the evolution Eq. (1) for the probability distribution of the meso- but regarding the large error bars introduced by our methods and
scopic stress: First (a) the yielding happens at a constant homoge- estimations it appears still convincing and suggest a very generic
neous yield stress σc . Then (b) the local response to surrounding flow behaviour.
plastic events is encompassed through local diffusive stress fluc- One of the strongest prediction of the HL model is that both
tuations with a well defined diffusion coefficient DHL (t). And the the prefactor A as well as the dynamical yield stress σy are de-
third assumption (c) concerns the relaxation of the local stress to termined solely by the specific value of the coupling constant be-
zero, once a site yields, leading to a typical stress jump ∆σHL ≈ σc . tween diffusion and plastic activity. In Fig. 3 we compare the
The HL model assumptions, aiming for a self-consistent descrip- dimensionless yield stress σy /hσcm i and prefactor A/(G0 σc τ)1/2 ,

1–8 | 5
 density α̃micro /10−1 α̃ f low /10−1
1.149 2.81 ± 0.21 3.21 ± 0.12
 

 

1.189 2.47 ± 0.27 3.30 ± 0.19
 
2.57 ± 0.26 3.28 ± 0.16
1.253 2.27 ± 0.24 

 3.28 ± 0.18 


1.322 2.74 ± 0.30 3.34 ± 0.18

Table 1 Comparison of different coupling constant measurements as explained in the text.

obtained from a Herschel Bulkley fit of the flow curves, with the original version of HL model equations such as the partial relax-
analytically obtained parametric curve. We observe that all the ation of the local stress after a yielding event and the introduction
data belonging to different density values collapse as expected of a yield stress distribution obtained from the microscopic simu-
by the previous collapse of the data roughly onto a single point lations. In a future work we plan to test the coherence of such a
close to the theoretical curve, that lies within the estimated er- modified model with our microscopic approach. Further it would
ror. This result points to a universal determination of the flow be highly desirable to test the degree of generality of the results
curves, a priori strongly sensitive to the density (note the density on other model systems 46 and experimental setups, that have ac-
dependence of A and σy in the inset of Fig. 3), through one single cess to the measure of local plastic activity 21,23–26 .
density independent parameter, namely the coupling constant α̃. Acknowledgments KM acknowledges financial support of the
As suggested by the KEP model 15 , we expect this parameter to be French Agence Nationale de la Recherche (ANR), under grant
only dependent on the specific form of the elastic propagator 9 . ANR-14-CE32-0005 (project FAPRES) and FP acknowledges fi-
If we assume the HL parametric relation between prefactor nancial support from ERC grant ADG20110209. Most of the
A(α̃) and dynamical yield stress σy (α̃) to hold, it is possible to computations were performed using the Froggy platform of the
estimate the coupling constant α̃ in an alternative way through CIMENT infrastructure supported by the Rhône-Alpes region
the macroscopic measure of the flow curve, in the following re- (GRANT CPER07-13 CIRA) and the Equip@Meso project (refer-
ferred to as α̃flow . Despite all the rough approximations we had to ence ANR-10-EQPX-29-01). Further we would like to thank Elis-
make, that tend to introduce large error bars on the data, we find abeth Agoritsas, Eric Bertin and Jean-Louis Barrat for fruitful dis-
the comparison between the coupling constant α̃micro with the al- cussions and valuable revisions of our manuscript. Also we would
ternative measurement using the expression for the rescaled yield like to thank Jörg Rottler for assistance in the cluster analysis of
stress† σy /σc (α̃) rather convincing (see table 1). Both measure- plastic events.
ments suggest a density independent result with a quite small
relative error of approximately 20%. Appendix A
Altogether our data suggests that revisiting the rule of setting
The parametric relation between the yield stress and the HB
the stress to zero after a yield event in the HL model equations, by
prefactor
changing the gain term in the evolution Eq. (1) together with the
introduction of a more realistic yield stress distribution, seems to In this section, we compute the relation between the dynamic
be a promising route to reach a more realistic mean-field model- yield stress σY and the prefactor A in the relation hσ i = σY + Aγ̇ 1/2
ing of athermally sheared amorphous systems. derived from the HL model (as given in Fig. 3). We use the
method developed in 38 : the probability density function solving
5 Conclusion the (stationary) HL equation (2) of the paper is expanded in the
In this study we aimed at testing some of the most basic assump- following way
tions and predictions of mean-field modeling for the rheology of
P(σ ) = Q0 (σ ) + γ̇ 1/2 Q1 (σ ) + · · · , for σ in [−σc , σc ] (9)
athermally sheared amorphous systems. In conclusion we obtain
a consistent picture of how to model correctly the mechanical
 
|σ | − σc
noise in the regime of large enough driving rates (far from the P(σ ) = γ̇ 1/2 R1± +··· , for ±σ in [σc , +∞[ (10)
γ̇ 1/2
true critical point 27 ). We find that we can incorporate the noise
into a normal diffusion of local stresses with a noise amplitude This ansatz has been proved to be correct in the case where
solely governed by the rate of plastic activity as proposed by the α < σc2 /2 (again see 38 ). Moreover, term-by-term integration is
Hébraud-Lequeux (HL) model 2 . We not only confirm this phys- allowed which means that
ical picture using molecular dynamics simulations, but we also
Z 
show that the coupling strength between diffusion and the rate of
Z
hσ i = σ P(σ )dσ = σ Q0 (σ )dσ
plastic activity, a dimensionless and density independent quantity,
| {z }
seems to determine the specific form of the rheological response. σY
Our data analysis suggest some important modifications in the Z 
+ σ Q1 (σ )dσ γ̇ 1/2 + · · · (11)
| {z }
† given in Appendix A A

6| 1–8
where the dots are terms of higher order than γ̇ 1/2 . Now all that In this system d˜ is an unknown coefficient to be simultaneously
is left to do is to identify Q0 and Q1 . By plugging (9) in (1) we computed with Q1 . The computation of Q1 is tedious but straight-
obtain the following equation on Q0 : forward. The expression of Q1 are quite lengthy but can be
 checked out on a symbolic computation program: if −σc ≤ σ ≤ 0
2 0 0 d
−d∂σ Q + τG0 ∂σ Q = α δ0 then


Q0 (±σc ) = 0 (12)
τG0 G0 τ
Q1 (σ ) = e d σ

R σc Q0 (σ )dσ = 1
 √  −G0 τσc
−σc σc d 1 − e d
Note that in this equation d is an unknown coefficient used to
˜ τ
dG G0 τ
enforce the integral condition on Q0 ; it is physically related to the −  0 −G0 τσc  (σ + σc ) e d σ
diffusion coefficient by τD ∼ d γ̇. d 2 σc 1 − e d
This system can easily be integrated and we find Q0 to be
1 (20)
+ G0 τ −G0 τ ×
e d σc −e d σc
 σ G τ −σ τG

1  e d0 − e cd 0 −σc ≤ σ ≤ 0
0 √ ! !
Q (σ ) = −σc G0 τ (σ −σc )G0 τ (13) σc G0 τ d˜ d d˜  −G0 τ 
σ
σc (1 − e d ) 1 − e d
 0 ≤ σ ≤ σc − + 1+e d c
×
d dσc σc2 dσc
and d is selected so that the following equation holds true:  G0 τ 
  e d (σ +σc ) − 1
α d σc G0 τ
= tanh (14)
σc2 σc G0 τ 2d and if 0 ≤ σ ≤ σc
Note that this equation (in d) has a unique solution d(α) if and τG0 G0 τ
Q1 (σ ) = √   e d (σ −σc )
only if 0 < α/σc2 < 1/2. However, in view of (13), it is easier −
G0 τσc
σc d 1 − e d
to express everything in terms of the parameter d instead of α:
then the limit α/σc2 → 1/2 is equivalent to d → +∞ and α → 0 is ˜ τ
dG G0 τ
equivalent to d → 0. +  0 G0 τσc  (σ − σc ) e d (σ −σc )
d 2 σc 1 − e− d
Now using the formula for computing σY we obtain by integra-
tion,   1
σY 1 σc G0 τ d + ×
= coth − (15) G0 τ G0 τ
σc 2 2d σc G0 τ e d σc − e− d σc

√ ! !
To obtain the prefactor A we must compute Q1 which necessi- σc G0 τ d˜ d d˜  G0 τ
σ

− − + 1+e d c
×
tates the computation of R1± . Using the continuity of ∂σ P at ±σc , d dσc σc2 dσc
one can find that R1± (z) are functions satisfying
 G0 τ 
 1 − e d (σ −σc )
−d∂z2 R1± + R1± = 0
(21)



∂ R1 (0) = ∂ Q0 (σ ) = −
 τG0
z + σ c −σc τG0
dσc (1−e d ) (16)

 −σc τG0 The parameter d˜ is selected to enforce the vanishing integral con-
−∂z R1− (0) = ∂σ Q0 (−σc ) = τG0 e σd c τG0

dition which amounts to taking

dσc (1−e− d )
τG0 σc τG0 σc
d 3/2 e d − e− d + d0 c
2τG σ
Again, it is easy to solve this system:
d˜ = − τG0 σc τG0 σc (22)
σc e d − e− d − 2τG0 σc
τG0 −z
√ d
R1+ (z) = √ −σc τG0 e
d (17)
dσc (1 − e d )
−σc τG0
τG0 e d −z
√
R1− (z) = √ σc τG0 e
d (18)
dσc (1 − e− d )
Finally we can compute the prefactor A which is equal to
σ Q1 (σ )dσ . All in all, we obtain (let us note u = (τG0 σc )/d):
R
Now we can write down the equations satisfied by Q1 using the
continuity of P at σ = ±σc : σY 1 1
= coth (u) − (23)
 σc 2 u
 ˜ 1 δ0 + ∂σ2 Q0 )
−d∂σ2 Q1 + G0 τ∂σ Q1 = d(
 α
1 eu (u cosh (2u) + 2 sinh (2u) + u)


Q 1 (σ ) = R1 (0) = τG0 A
√ = √


 c + √ −σc τG0 (24)
2 u u (sinh (2u) − u)

dσc (1−e
d ) σc τG0
−σc τG0 (19)
τG0 e d
Q1 (−σc ) = R1− (0) = √ 6u2 coth (2u)
 σc τG0

 dσc (1−e− d ) 1
− √

R
2 u u (sinh (2u) − u)

 σc 1
−σc Q (σ )dσ = 0

1–8 | 7
References 25 A. Le Bouil, A. Amon, S. McNamara and J. Crassous, Physical
Review Letters, 2014, 112, 246001.
1 P. Sollich, F. Lequeux, P. Hebraud and M. E. Cates, Physical
26 K. E. Jensen, D. A. Weitz and F. Spaepen, Physical Review E,
Review Letters, 1997, 78, 2020–2023.
2014, 90, 042305.
2 P. Hébraud and F. Lequeux, Physical Review Letters, 1998, 81,
27 C. Liu, E. Ferrero, F. Puosi, J. Barrat and K. Martens,
2934–2937.
arXiv:1506.08161, 2015.
3 M. L. Falk and J. S. Langer, Phys. Rev. E, 1998, 57, 7192–7205.
28 M. Kac, G. E. Uhlenbeck and P. C. Hemmer, J. Math. Phys.,
4 D. Rodney, A. Tanguy and D. Vandembroucq, Modelling Simul.
1963, 4, 216.
Mater. Sci. Eng., 2011, 19, 083001. 29 W. Klein, H. Gould, N. Gulbahce, J. B. Rundle and K. Tiampo,
5 A. Nicolas, K. Martens and J. L. Barrat, Epl, 2014, 107, 6. Physical Review E, 2007, 75, 031114.
6 E. Agoritsas, E. Bertin, K. Martens and J.-L. Barrat, The Euro- 30 N. Xu and C. O’Hern, Physical Review Letters, 2014, 94,
pean Physical Journal E, 2015, 38, 71. 055701.
7 A. S. Argon, Acta Metallurgica, 1979, 27, 47–58. 31 J. Puckett and K. Daniels, Physical Review Letters, 2013, 110,
8 J. D. Eshelby, Proc. R. Soc. London A, 1957, 241, 376. 058001.
9 F. Puosi, J. Rottler and J. L. Barrat, Physical Review E, 2014, 32 J. Langer, arXiv:1501.07228, 2015.
89, 042302. 33 S. Karmakar, E. Lerner, I. Procaccia and J. Zylberg, Physical
10 A. Lemaître and C. Caroli, Phys. Rev. Lett., 2009, 103, 065501. Review E, 2010, 82, 031301.
11 K. Martens, L. Bocquet and J.-L. Barrat, Phys. Rev. Lett., 2011, 34 V. Mansard, A. Colin, P. Chauduri and L. Bocquet, Soft Matter,
106, 156001. 2011, 7, 5524–5527.
12 K. M. Salerno and M. O. Robbins, Physical Review E, 2013, 88, 35 Y. Gati, PhD thesis, Ecole des ponts ParisTech, 2004.
062206. 36 E. Cances, I. Catto and Y. Gati, SIAM J. Math. Anal., 2006, 37,
13 Z. Budrikis and S. Zapperi, Physical Review E, 2013, 88, 60.
062403. 37 E. Cances, I. Catto, Y. Gati and C. L. Bris, Multiscalce Model.
14 J. Lin, E. Lerner, A. Rosso and M. Wyart, Proceedings of the Simul., 2006, 4, 1041.
National Academy of Sciences of the United States of America, 38 J. Olivier, PhD thesis, Université de Grenoble, 2011.
2014, 111, 14382–14387. 39 J. Olivier, Z. Angew. Math. Phys., 2010, 61, 445.
15 L. Bocquet, A. Colin and A. Ajdari, Physical Review Letters, 40 J. Olivier and M. Renardy, SIAM J. Appl. Math., 2011, 71,
2009, 103, 036001. 1144.
16 K. A. Dahmen, Y. Ben-Zion and J. T. Uhl, Nature Physics, 2011, 41 J. Olivier, Science China-Mathematics, 2012, 55, 435–452.
7, 554–557.
42 R. Brüning, D. A. St-Onge, S. Patterson and W. Kob, J. Phys.:
17 I. Cantat and O. Pitois, Journal of Physics: Condensed Matter, Condens. Matter, 2009, 21, 035117.
2005, 17, S3455.
43 T. Schneider and E. Stoll, Phys. Rev. B, 1978, 17, 1302–1322.
18 L. Bécu, S. Manneville and A. Colin, Physical review letters,
44 K. Salerno, C. Maloney and M. Robbins, Physical Review Let-
2006, 96, 138302.
ters, 2012, 109, 105703.
19 J. Goyon, A. Colin, G. Overlez, A. Ajdari and L. Bocquet, Na-
45 A. Nicolas, J. Rottler and J.-L. Barrat, Eur. Phys. J. E Soft Mat-
ture, 2008, 454, 84–87.
ter, 2014, 37, 37–50.
20 J. Goyon, A. Colin and L. Bocquet, Soft Matter, 2010, 6, 2668–
46 C. Fusco, T. Albaret and A. Tanguy, European Physical Journal
2678.
E, 2014, 37, 43.
21 A. Amon, V. B. Nguyen, A. Bruand, J. Crassous and
47 A. Tanguy, F. Leonforte and J.-L. Barrat, Eur. Phys. J. E, 2006,
E. Clement, Physical Review Letters, 2012, 108, 135502.
20, 355–364.
22 K. Kamrin and G. Koval, Physical Review Letters, 2012, 108,
48 P. Sollich and A. Barra, In preparation.
178301.
49 P. Sollich, School on Glass Formers and Glasses, Bangalore,
23 K. W. Desmond, P. J. Young, D. D. Chen and E. R. Weeks, Soft
India, 2010.
Matter, 2013, 9, 3424–3436.
50 H. Mizuno, S. Mossa and J.-L. Barrat, Phys. Rev. E, 2013, 87,
24 E. D. Knowlton, D. J. Pine and L. Cipelletti, Soft Matter, 2014,
042306.
36, 6931–6940.

8| 1–8