<figure><img src="image_1.png" /></figure>
<figure><img src="image_1.png" /></figure>

<figure><img src="image_2.png" /></figure>
arXiv
<figure><img src="image_1.png" /></figure>
<figure><img src="image_1.png" /></figure>

<figure><img src="image_2.png" /></figure>
# Probing relevant ingredients in mean-field approaches for the athermal rheology of yield stress materials

Francesco Puosi, Julien Olivier, Kirsten Martens

► To cite this version:

Francesco Puosi, Julien Olivier, Kirsten Martens. Probing relevant ingredients in mean-field approaches for the athermal rheology of yield stress materials. Soft Matter, 2015, 11 (38), pp.7639-7647. 10.1039/C5SM01694K . hal-01263149

HAL Id: hal-01263149

https://hal.science/hal-01263149

Submitted on 1 Feb 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Regular paper
<figure><img src="image_1.png" /></figure>

<figure><img src="image_2.png" /></figure>
<figure><img src="image_2.png" /></figure>

CrossMark

<figure><img src="image_3.png" /></figure>

Cite this: *Soft Matter*, 2015, 11, 7639

Received 9th July 2015,
Accepted 3rd August 2015

DOI: 10.1039/c5sm01694k

www.rsc.org/softmatter

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

Francesco Puosi,*abc Julien Olivier^d and Kirsten Martens^bc

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.

## 1 Introduction

The theoretical understanding of the yielding transition in athermally driven disordered systems is a highly challenging problem and no consensus has been established even on the basic ingredients that should underlie coarse grained descriptions of the non-linear rheological response of yield stress materials.<sup>1–6</sup>

The only commonly accepted point of view is that disordered materials, such as glasses or soft matter, exhibit a strongly heterogeneous 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 elastically.<sup>7</sup> These plastic events induce long-range elastic deformations in the system leading to complex correlations of the yielding regions in form of plastic avalanches.<sup>8–14</sup>

One of the major concerns remains to know whether despite this dynamical complexity it is sensible to describe the yielding dynamics within mean-field descriptions.<sup>1–3,15,16</sup> We tackle this key question using particle based simulations, concentrating on flow responses in a regime that is relevant for many rheological setups.<sup>17–27</sup> We find that at high enough shear rates and/or small enough system sizes we recover a dynamics

well described by mean-field considerations,<sup>28</sup> similar to near mean-field critical points in equilibrium phase transitions.<sup>29,30</sup>

An important feature that found its way into several theories, e.g. the Soft Glassy Rheology (SGR) and the Shear Transformation Zone Theory (STZ), is to describe yielding in a generalized thermodynamic description with an effective temperature.<sup>1,3</sup> But the analogy with thermally activated events proposed by the SGR theory has recently been questioned for athermal rheology.<sup>5,6</sup> Also it has been shown that the thermodynamic interpretation of the effective temperature is problematic in the athermal regime<sup>31,32</sup> and that the STZ theory is not able to predict the non-linear response,<sup>33</sup> expected in the small driving limit.<sup>28</sup>

In this work we show that the dynamically created noise in sheared amorphous systems can be encompassed through a normal diffusion in the local stresses<sup>34</sup> (for a schematic view see Fig. 1).

<figure><img src="image_4.png" /><figcaption>**Fig. 1** Schematic view of the coarse grained picture: the particle simulation box is divided into smaller mesoscopic parts that can hold exactly one plastic event. The plastic activity in the surrounding of a mesoscopic region (see graphical interpretation on the right) leads to stress fluctuations around the mean value of the stress (dashed line in the graph on the right), which is controlled by the external forcing.</figcaption></figure>

<footnote><sup>a</sup> 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.
E-mail: francesco.puosi@ens-lyon.fr</footnote>

<footnote><sup>b</sup> Université Grenoble Alpes, LIPHY, F-38000 Grenoble, France</footnote>

<footnote><sup>c</sup> CNRS, LIPHY, F-38000 Grenoble, France</footnote>

<footnote><sup>d</sup> Aix Marseille Université, Centre Mathématiques et Informatique, Technôpole Château Gombert, 39 Rue F Joliot Curie, F-13453 Marseille 13, France</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 2 Flow curves: rescaled macroscopic shear stress $\sigma_{xy}/\sigma_c$ as a function of the rescaled applied shear rate $\dot{\gamma}(G_0\tau_0/\sigma_c)$ for different densities (see text). The dashed line is a fit with the Herschel-Bulkley form $\sigma_{xy} = \sigma_\gamma + A\dot{\gamma}^{1/2}$. The inset shows the same data before rescaling, dashed lines display a guide to the eye with the same Herschel-Bulkley fitting form.</figcaption></figure>

In agreement with former works,⁵,⁶ we consider that local stress diffusion is acting as noise in the tilt of the local potential energy landscape with an amplitude proportional to the plastic activity, an observable that can be measured experimentally.²¹,²⁴

We analyse within molecular dynamics simulations the coupling strength between the mechanical noise and the plastic activity. Within our simulations we successfully relate this coupling strength, a dimensionless and density independent quantity, in a consistent manner to the flow response without any further parameters. This finding leads us to reconsider a speculative mean-field scenario put forward a long time ago in the so-called Hébraud-Lequeux (HL) model,² that was at the basis of important further developments,¹⁵,²²,³⁵ providing so far one of the best self-consistent mean-field description of mechanical noise in athermally sheared disordered systems.

## 2 Mean-field approach

Within the framework of the HL-model we can establish a link between the dynamical yield stress and the prefactor of a Herschel Bulkley type power law fit of the rheological curves of athermally sheared yield stress materials. This relation should, according to this theory, solely be determined by the coupling strength between mechanical noise and the rate of activity. In the following we review briefly the main assumptions and results of the mean-field approach.² The model describes the state of a soft glassy material via the probability density $\mathcal{P}$ of local shear stresses $\sigma = \sigma_{xy}$ in regions of mesoscopic size W while the material is sheared at rate $\dot{\gamma}$. The time evolution of $\mathcal{P}$ is given by

$$ \partial_t \mathcal{P}(\sigma, t) = -G_0 \dot{\gamma}(t) \partial_\sigma \mathcal{P} - \frac{1}{\tau} \theta(|\sigma| - \sigma_c) \mathcal{P} \\ + \Gamma(t) \delta(\sigma) + D_{\text{HL}}(t) \partial_\sigma^2 \mathcal{P} \quad (1) $$

where $\theta(x)$ and $\delta(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 $\partial_\sigma \mathcal{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 $\sigma_c$ at a rate given by $1/\tau$. It has been argued that the phenomenological parameter $\tau$ can be interpreted as the duration of a plastic event in the low shear rate limit.⁶ 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

$$ \Gamma(t) = \frac{1}{\tau} \int_{|\sigma| > \sigma_c} \mathcal{P}(\sigma, t) d\sigma. \quad (2) $$

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 dimensionless coupling constant $\tilde{\alpha}$

$$ D_{\text{HL}}(t) = \tilde{\alpha} \sigma_c^2 \Gamma(t). \quad (3) $$

This last relation introduces a non-linearity into eqn (1), since the rate of plastic activity itself depends on the density probability of the mesoscopic stresses. It is this coupling that renders the problem non-trivial and yields interesting results regarding the behaviour of macroscopic quantities.

This model is known to exhibit a unique stationary state for a finite shear rate in the large time limit, where the probability density for the mesoscopic stresses becomes time independent. To determine the time averaged macroscopic stress in the steady state one averages over the mesoscopic stresses weighted by the corresponding steady state probability density

$$ \langle \sigma \rangle = \int \sigma \mathcal{P}(\sigma) d\sigma. \quad (4) $$

Using appropriate units we can write the equations in dimensionless quantities, expressing stress related values in units of the local yield stress $\sigma_c$, time quantities in units of $\tau$, the shear rate in units of $\sigma_c/(G_0\tau)$ and the stress diffusion coefficient in units of $\sigma_c^2/\tau$, leaving only two independent dimensionless model parameters that determine the flow behaviour, namely the dimensionless shear rate and coupling constant $\tilde{\alpha}$.

The rheological results in the small shear rate limit for this model are well studied.²,⁶,³⁶⁻⁴² For small enough coupling strength $\tilde{\alpha} < 1/2$ the HL model predicts a Herschel-Bulkley flow behaviour of exponent $1/2$, $\langle \sigma \rangle \approx \sigma_Y + A\dot{\gamma}^{1/2}$, for the time averaged macroscopic stress in the steady state $\langle \sigma \rangle = \int \sigma \mathcal{P}(\sigma) d\sigma$, with the two constants $\sigma_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 $\tilde{\alpha}$. This means that we will obtain a universal curve† for the rescaled quantities $A/\sqrt{G_0\tau\sigma_c}$ versus $\sigma_Y/\sigma_c$ parametrized through $\tilde{\alpha}$ (see Fig. 3).

<footnote>† For the analytical derivation of this constant relation see the Appendix A</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 3 HL predictions vs. MD simulations: (a) the points are the rescaled Herschel Bulkley prefactor A as a function of the rescaled dynamical yield stress σY for four different densities, error bars are estimated from the measurements of the different quantities. The full line corresponds 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 from a very simplified picture.</figcaption></figure>

The physical interpretation of this coupling constant could be a mechanical fragility² that would depend on the details of the microscopic interactions and thus should be material dependent. In a later work by Bocquet et al., the authors use a spatial approach to the problem that allows for the derivation of an expression for the coupling constant $\tilde{\alpha}$ if one assumes a decorrelation of the plastic event dynamics. Within this approximation it is possible to express $\tilde{\alpha}$ as a function of the specific form of the elastic stress propagator.¹⁵ This suggests that there should be classes of materials with similar elastic responses that not only share the same non-linear exponent in the flow curve, but also comparable relations between the dynamical yield stress and the prefactor in the Herschel Bulkley form.

In the following we aim at probing not only the above constant relation using microscopic dynamic simulations, but also to test 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.

Quantities and parameters like the local stress and the shear modulus $G_0$ are well defined quantities in the microscopic simulations and rather easy to measure. Also, if stress fluctuations turn out to be diffusive, $D_{HL}$ has a well defined meaning. However, other HL parameters like the local yield stress $\sigma_c$, the phenomenological parameter $\tau$ and the rate of plastic activity $\Gamma(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.

## 3 Microscopic model

We have investigated a generic two-dimensional (2D) model of a 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_{\alpha\beta}(r) = 4\epsilon_{\alpha\beta}[(\sigma_{\alpha\beta}/r)^{12} - (\sigma_{\alpha\beta}/r)^6]$ with $\alpha, \beta = A, B$ and $r$ being the distance between two particles. The parameters $\epsilon_{AA}$, $\sigma_{AA}$ and $m_A$ define the units of energy, length and mass; the unit of time is given by $\tau_0 = \sigma_{AA}\sqrt{(m_A/\epsilon_{AA})}$. We set $\epsilon_{AA} = 1.0$, $\epsilon_{AB} = 1.5$, $\epsilon_{BB} = 0.5$, $\sigma_{AA} = 1.0$, $\sigma_{AB} = 0.8$ and $\sigma_{BB} = 0.88$ and $m_A = m_B = 1$. This choice is known to prevent crystallization in 2D at low temperature.⁴³ 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 $\delta t = 0.005$. The athermal limit is achieved by thermostating the system at zero temperature via a Langevin thermostat⁴⁴ with a damping coefficient $\zeta = 1$ which corresponds to a strongly overdamped condition for the dynamics.⁴⁵ This model, which has been widely studied in different versions,³,¹⁰,⁴⁵,⁴⁶ is usually considered appropriate for colloidal and other soft glasses. This choice is motivated by the purpose of investigating the general aspects of the rheology of athermal yield stress materials.

To investigate different athermal flow responses, we explore states with different number density $\rho = (N_A + N_B)/V$ by changing the volume $V$ of the system. 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 rate $\dot{\gamma}$ by deforming the box dimensions and remapping the particle positions. The quenching protocol has virtually no effect since we focus on the steady state shear (total imposed deformation $\Delta\dot{\gamma} > 20\%$).

In order to characterize the plastic activity of the system we consider the $D_{min}^2$ quantity.³ For a given particle $i$, $D_{min}^2$ is defined as the minimum over all possible linear deformation tensors $\epsilon$ of:

$$D^2(i, t, \delta t) = \sum_j [\mathbf{r}_{ij}(t + \delta t) - (\mathbf{I} + \epsilon) \cdot \mathbf{r}_{ij}(t)]^2 \quad (5)$$

where the index $j$ runs over all the neighbors of the reference particle $i$ and $\mathbf{I}$ is the identity matrix. We set the time lag to $\delta t = 4$. This value is a compromise between having a good signal, i.e. large irreversible displacements, and being able to resolve individual plastic events.

### 3.1 Macroscopic flow curve

In Fig. 2 we show the dependence of the macroscopic shear stress $\sigma_{xy}$ on the applied shear rate $\dot{\gamma}$. The flow curves are well described by the Herschel-Bulkley (HB) form, $\sigma_{xy} = \sigma_Y + A\dot{\gamma}^n$, with an exponent $n \approx 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 $\tilde{\alpha} < 1/2$, gives indistinguishable results. Although other works on sheared disordered material in two and three dimensions report similar values for the flow curve exponent,⁴⁶,⁴⁷ 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.¹³,¹⁴,²⁸,⁴⁵ 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.

## 3.2 Size of an elementary plastic region

Here we describe the procedure we followed to convert the microscopic simulations into a mesoscopic description. First, we denote a particle i as active, i.e. performing a plastic rearrangement, at a given time t if the corresponding $D_{min}^2(i,t)$ is larger than a threshold value that we fix equal to 0.1. In Fig. 4(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 is reasonable. Then a high-resolution discretization of the system is performed by dividing it into $k \times k$ square blocks of length $w_0 = 2$. A small block is considered as active if it contains at least one active particle.

The mean size, i.e. the mean linear extension, of a plastic event $\langle l_{pl} \rangle$ can be estimated by a cluster analysis of the spatial arrangement of the active blocks in the configurations explored by the system. We employed a modified version of the Hoshen-Kopelman algorithm$^{48}$ in order to account for periodic boundary conditions. If we assume $\langle l_{pl} \rangle = \langle A_{pl} \rangle^{1/2}$ where $A_{pl}$ is the area of a plastic cluster, we obtain $\langle l_{pl} \rangle \approx 6$, with no relevant dependence on the density and shear rate (in the range $\dot{\gamma} \le 10^{-4}$). The value is in accord with previous works reporting plastic regions with a size of a few particle diameters.$^{46,49}$ Furthermore, this agreement justifies the criterion we employed to define active particles. Indeed, if we improve the resolution in the analysis of the plastic

<figure><img src="image_7.png" /><figcaption>Fig. 4 Activity and stress diffusion: (a) the full line represent a part of the macroscopic stress-strain curve in the steady state regime for a density $\rho = 1.149$ and a shear rate $\dot{\gamma} = 10^{-6}$. The average slope of the elastic parts on the curve yields $G_0$. Open symbols: corresponding evolution of the number of active particles $N_{pl}$ in the system. (b) Coarse mean-square stress difference $\langle (\Delta\sigma_{xy})^2(t) \rangle$ as a function of time for different values of the density. Note the linear behaviour at short times. Inset: Stress diffusion coefficient $D_\sigma = \langle (\Delta\sigma_{xy})^2(t) \rangle/t$ as a function of density.</figcaption></figure>

activity by decreasing the threshold on $D_{min}^2$ by a factor 10, while the number of active particles increase by a factor 4, the mean extension of a plastic event $\langle l_{pl} \rangle$ is reduced by half, suggesting that single particle rearrangements are erroneously taken into account.

Next we implement the coarse-graining of microscopic simulations on the scale of individual plastic events. The simulation box is divided into $M \times M$ square blocks with $M$ chosen in order to have $W \approx \langle l_{pl} \rangle$. The local shear stress $\sigma_{xy}^m$ of a block m is defined as:

$$ \sigma_{xy}^m = -\frac{1}{W^2} \sum_{i \in m} m v_{i,x} v_{i,y} + \frac{1}{2W^2} \sum_{i \in m} \sum_{j=1}^{N} \frac{\partial V(r^{ij}) r_{x,y}^{ij}}{\partial r^{ij}} \quad (6) $$

where $v_{i,x}$ and $v_{i,y}$ are the x and y components of the velocity of the particle i, $r^{ij}$ is the distance between the particles i and j and the summation of i is performed over the particles in the block. The macroscopic stress tensor $\sigma_{xy}$ is obtained by the summation of $\sigma_{xy}^m$ over all the blocks.

## 3.3 Stress diffusion and duration of a plastic event

The above introduced mean-field model assumes local stress fluctuations, obeying a normal diffusion process in stress space. To test this idea we define the coarse mean-square stress difference as:

$$ \langle (\Delta\sigma_{xy}^m)^2(t) \rangle = \langle (\sigma_{xy}^m(t_0 + t) - \sigma_{xy}^m(t_0) - G_0\dot{\gamma}t)^2 \rangle \quad (7) $$

where $\sigma_{xy}^m(t)$ is the stress in a given block at time t and the last term in eqn (7) accounts for the stress increase due to the elastic deformation of the system, being $G_0$ the macroscopic shear modulus. In Fig. 4(b), the short time behavior of $\langle (\Delta\sigma_{xy}^m)^2(t) \rangle$ is shown for a finite shear rate $\dot{\gamma} = 10^{-6}$, approaching the quasistatic limit. We observe that at short times, $\langle (\Delta\sigma_{xy}^m)^2(t) \rangle$ increases linearly with time. We define the stress diffusion coefficient as $D_\sigma = \langle (\Delta\sigma_{xy}^m)^2(t) \rangle/t$. In the inset of Fig. 4(b) we show the dependence of $D_\sigma$ on the density of the system.

To define a duration of a plastic event we analyse the two-time autocorrelation function of the $D_{min}^2$ quantity for active particles. In Fig. 5(a) we show this correlation function $C_p = \langle D_{min}^2(t_0) D_{min}^2(t_0 + t) \rangle / \langle (D_{min}^2(t_0))^2 \rangle$ as a function of time. We observe that $C_p$ decays exponentially with a characteristic time $\tau_p$, that depends weakly on density (see the inset of Fig. 5(a)) and that is close to the damping time $\tau_d = \xi^{-1}$ of the Langevin thermostat. We choose to interpret this decorrelation time as the typical duration of a plastic event entering the HL model description.

## 3.4 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.$^{50–52}$ The system is frozen except for a target region, i.e., a mesoscopic block, and it's subjected to a simple shear deformation, with a quasi-static protocol. The frozen region can only deform affinely whereas
<figure><img src="image_8.png" /><figcaption>Fig. 5 Duration of events and local yield stress: (a) $D_{\text{min}}^2$ autocorrelation function for active particles $C_p = (D_{\text{min}}^2(t_0)D_{\text{min}}^2(t_0+t))/(D_{\text{min}}^2(t_0))^2$ as a function of time. Lines are fitting curves $\exp(-t/\tau_p)$. Inset: Dependence of $\tau_p$ on the density. (b) Distribution of the local yield stress $p(\sigma_c^m)$ for different values of the density. Full lines are Gaussian distributions fits. Left inset: Mean value $\langle\sigma_c^m\rangle$ and variance $\langle\delta\sigma_c^m\rangle$ as a function of density. Dashed lines are power-law guides. Right inset: Absolute local stress drop $\langle\Delta\sigma^m\rangle$ and the relative local stress drop $\langle\Delta\sigma^m\rangle/(\sigma_c^m)$ as a function of density. Dashed lines are power-law and constant value guides to the eye. Symbol code is the same as in panel (a).</figcaption></figure>

the target block is allowed to relax non-affinely. For small deformations the target region behaves elastically and the stress increases linearly, with a slope controlled by the local shear modulus.⁵² 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 stress-strain curve. For a given block, we define as the local yield stress $\sigma_c^m$ the value of the local stress at the first maximum.

In Fig. 5(b) we show the distribution of $\sigma_c^m$ 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 $\sigma_c$. The exposed distributions are well described by a simple Gaussian forms. The mean values $\langle\sigma_c^m\rangle$ and the variance $\langle\delta\sigma_c^m\rangle$ of the distributions are shown in the inset of the figure. As the density increases, $\langle\sigma_c^m\rangle$ increases, due to the enhancement of the repulsive interactions between particles.

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 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 $\langle\Delta\sigma^m\rangle/\langle\sigma_c^m\rangle \approx 0.2$ (see inset of Fig. 5).

In a former work⁵² 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 \approx 6$ the discrepancy is up to 50%, a significant error. We are aware that the estimate of $\langle\sigma_c^m\rangle$ 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,⁶ but this is beyond the scope of this study and left for further investigations.

## 4 Robustness of the HL model

In the following we aim at probing the robustness of the HL model by comparing our data from the microscopic simulations to the various assumptions and predictions of the mean-field description.

Let us first recall some of the basic assumptions made to write the evolution eqn (1) for the probability distribution of the mesoscopic stress: first (a) the yielding happens at a constant homogeneous yield stress $\sigma_c$. Then (b) the local response to surrounding plastic events is encompassed through local diffusive stress fluctuations with a well defined diffusion coefficient $D_{\text{HL}}(t)$. And the third assumption (c) concerns the relaxation of the local stress to zero, once a site yields, leading to a typical stress jump $\Delta\sigma_{\text{HL}} \approx \sigma_c$. The HL model assumptions, aiming for a self-consistent description of the mechanical noise, thus predict a diffusion coefficient proportional to the square of this typical local stress jump $\Delta\sigma_{\text{HL}}$ divided by a typical time scale, given by the inverse plasticity rate $\Gamma(t)$.

(a) Our data analysis reveals that a homogeneous yield stress is of course not verified in a disordered system, where one expects a distribution of yield stresses (see Fig. 5b). However, it has been recently shown that the existence of a local yield stress distribution in the HL dynamics does not strongly alter the predictions for the flow behaviour.⁶

(b) We tested as well the second assumption of a normal diffusion of the mesoscopic stresses, and we find that within the range of shear rates that we consider, the measurement of the average mean-square stress differences indeed allows for the determination of a well-defined diffusion coefficient. The HL coupling constant $\tilde{\alpha}_{\text{micro}}^{\text{HL}} = D_{\sigma}/(\sigma_c^2\Gamma)$ (see eqn (3)) with $D_{\sigma}$ from the results in Fig. 4(b) and $\Gamma$ from the cluster analysis turns out to be very small, with $\tilde{\alpha}_{\text{micro}}^{\text{HL}}$ of order $10^{-2}$ for all considered densities.

(c) However, our study suggests that the partial relaxation of the local stresses after yielding, $\langle\Delta\sigma^m\rangle/\langle\sigma_c^m\rangle \approx 0.2$ (instead of the assumed total relaxation) introduces an important corrective factor for the coupling constant. Since the typical stress jump in the diffusion process is now reduced, we obtain for the stress diffusion coefficient an altered expression

$$D_{\sigma} = \tilde{\alpha}_{\text{micro}}(\Delta\sigma^{m})^{2}\Gamma. \qquad (8)$$
<table><caption>**Table 1** Comparison of different coupling constant measurements as explained in the text</caption><thead><tr><th>Density</th><th>α̃<sub>micro</sub>/10<sup>-1</sup></th><th></th><th>α̃<sub>flow</sub>/10<sup>-1</sup></th><th></th></tr></thead><tbody><tr><td>1.149</td><td>2.81 ± 0.21</td><td rowspan="4">2.57 ± 0.26</td><td>3.21 ± 0.12</td><td rowspan="4">3.28 ± 0.16</td></tr><tr><td>1.189</td><td>2.47 ± 0.27</td><td>3.30 ± 0.19</td></tr><tr><td>1.253</td><td>2.27 ± 0.24</td><td>3.28 ± 0.18</td></tr><tr><td>1.322</td><td>2.74 ± 0.30</td><td>3.34 ± 0.18</td></tr></tbody></table>

We report the values obtained for this new definition of the coupling constant in Table 1. It turns out that these values are in good approximation density independent and yield an average value of about $\tilde{\alpha} \approx 0.26$. Thus the coupling constant displays a reasonable value smaller than 0.5 for which the HL model predicts⁶ a well-established Herschel Bulkley regime in the flow curve with an exponent of 1/2.

And indeed our flow curves from the microscopic simulations can be well fitted with an Herschel Bulkley expression $\sigma_{xy} = \sigma_y + A\tilde{\gamma}^{1/2}$ (see inset of Fig. 2). In Fig. 2 we plot our data in dimensionless units, rescaled as suggested by the HL model, using dimensionless quantities $\sigma_{xy}/\sigma_c$ and $\tilde{\gamma}(G_0\tau/\sigma_c)$. The shear modulus $G_0$ can be easily accessed in simulations as the slope of the stress–strain curve in the elastic regime. We approximate the critical stress $\sigma_c$ as the mean value $\langle\sigma_c^m\rangle$ of the local yield stress obtained with the frozen matrix method and we associate the phenomenological parameter $\tau$ with the duration of a plastic event in the low shear rate limit.⁶ Hence, we will approximate $\tau$ by the value of $\tau_p$ obtained through the measurement of the $D_{\min}^2$ two-time autocorrelation function. We observe a collapse of the flow curves for different densities. This collapse is not perfect, but regarding the large error bars introduced by our methods and estimations it appears still convincing and suggest a very generic flow behaviour.

One of the strongest prediction of the HL model is that both the prefactor $A$ as well as the dynamical yield stress $\sigma_y$ are determined solely by the specific value of the coupling constant between diffusion and plastic activity. In Fig. 3 we compare the dimensionless yield stress $\sigma_y/\langle\sigma_c^m\rangle$ and prefactor $A/(G_0\sigma_c\tilde{\gamma})^{1/2}$, obtained from a Herschel Bulkley fit of the flow curves, with the analytically obtained parametric curve. We observe that all the data belonging to different density values collapse as expected by the previous collapse of the data roughly onto a single point close to the theoretical curve, that lies within the estimated error. This result points to a universal determination of the flow curves, a priori strongly sensitive to the density (note the density dependence of $A$ and $\sigma_y$ in the inset of Fig. 3), through one single density independent parameter, namely the coupling constant $\tilde{\alpha}$. As suggested by the KEP model,¹⁵ we expect this parameter to be only dependent on the specific form of the elastic propagator.⁹

If we assume the HL parametric relation between prefactor $A(\tilde{\alpha})$ and dynamical yield stress $\sigma_y(\tilde{\alpha})$ to hold, it is possible to estimate the coupling constant $\tilde{\alpha}$ in an alternative way through the macroscopic measure of the flow curve, in the following referred to as $\tilde{\alpha}_{flow}$. Despite all the rough approximations

we had to make, that tend to introduce large error bars on the data, we find the comparison between the coupling constant $\tilde{\alpha}_{micro}$ with the alternative measurement using the expression for the rescaled yield stress$‡$ $\sigma_y/\sigma_c(\tilde{\alpha})$ rather convincing (see Table 1). Both measurements suggest a density independent result with a quite small relative error of approximately 20%.

Altogether our data suggests that revisiting the rule of setting the stress to zero after a yield event in the HL model equations, by changing the gain term in the evolution eqn (1) together with the introduction of a more realistic yield stress distribution, seems to be a promising route to reach a more realistic mean-field modeling of athermally sheared amorphous systems.

## 5 Conclusion

In this study we aimed at testing some of the most basic assumptions and predictions of mean-field modeling for the rheology of athermally sheared amorphous systems. In conclusion we obtain a consistent picture of how to model correctly the mechanical noise in the regime of large enough driving rates (far from the true critical point$^{28}$). We find that we can incorporate the noise into a normal diffusion of local stresses with a noise amplitude solely governed by the rate of plastic activity as proposed by the Hébraud–Lequeux (HL) model.² We not only confirm this physical picture using molecular dynamics simulations, but we also show that the coupling strength between diffusion and the rate of plastic activity, a dimensionless and density independent quantity, seems to determine the specific form of the rheological response. Our data analysis suggest some important modifications in the original version of HL model equations such as the partial relaxation of the local stress after a yielding event and the introduction of a yield stress distribution obtained from the microscopic simulations. In a future work we plan to test the coherence of such a modified model with our microscopic approach. Further it would be highly desirable to test the degree of generality of the results on other model systems⁴⁷ and experimental setups, that have access to the measure of local plastic activity.²¹,²³⁻²⁶

<footnote>‡ Given in Appendix A</footnote>
# Appendix A

## The parametric relation between the yield stress and the HB prefactor

In this section, we compute the relation between the dynamic yield stress $\sigma_Y$ and the prefactor $A$ in the relation $\langle\sigma\rangle = \sigma_Y + A\dot{\gamma}^{1/2}$ derived from the HL model (as given in Fig. 3). We use the method developed in:³⁹ the probability density function solving the (stationary) HL eqn (2) of the paper is expanded in the following way

$$ \mathcal{P}(\sigma) = Q^0(\sigma) + \dot{\gamma}^{1/2}Q^1(\sigma) + \dots, \quad \text{for } \sigma \text{ in } [-\sigma_c, \sigma_c] \quad (9) $$

$$ \mathcal{P}(\sigma) = \dot{\gamma}^{1/2} R_{\pm}^{1} \left( \frac{|\sigma| - \sigma_c}{\dot{\gamma}^{1/2}} \right) + \dots, \quad \text{for } \pm \sigma \text{ in } [\sigma_c, +\infty) \quad (10) $$

This ansatz has been proved to be correct in the case where $\alpha < \sigma_c^2/2$ (again see ref. 39). Moreover, term-by-term integration is allowed which means that

$$ \begin{aligned} \langle \sigma \rangle &= \int \sigma \mathcal{P}(\sigma) d\sigma \\ &= \underbrace{\left( \int \sigma Q^0(\sigma) d\sigma \right)}_{\sigma_Y} + \underbrace{\left( \int \sigma Q^1(\sigma) d\sigma \right)}_{A} \dot{\gamma}^{1/2} + \dots \end{aligned} \quad (11) $$

where the dots are terms of higher order than $\dot{\gamma}^{1/2}$. Now all that is left to do is to identify $Q^0$ and $Q^1$. By plugging (9) in (1) we obtain the following equation on $Q^0$:

$$ \begin{cases} -d\partial_{\sigma^2}Q^0 + \tau G_0 \partial_{\sigma}Q^0 = \frac{d}{\alpha}\delta_0 \\ Q^0(\pm\sigma_c) = 0 \\ \int_{-\sigma_c}^{\sigma_c} Q^0(\sigma)d\sigma = 1 \end{cases} \quad (12) $$

Note that in this equation $d$ is an unknown coefficient used to enforce the integral condition on $Q^0$; it is physically related to the diffusion coefficient by $\tau D \sim d\dot{\gamma}$.

This system can easily be integrated and we find $Q^0$ to be

$$ Q^0(\sigma) = \frac{1}{\sigma_c \left( 1 - e^{-\frac{\sigma_c G_0 \tau}{d}} \right)} \begin{cases} \left( \frac{\sigma G_0 \tau}{e^d} - \frac{-\sigma_c \tau G_0}{e^d} \right) - \sigma_c \le \sigma \le 0 \\ \frac{(\sigma - \sigma_c) G_0 \tau}{1 - e^d} \quad 0 \le \sigma \le \sigma_c \\ \end{cases} \quad (13) $$

and $d$ is selected so that the following equation holds true:

$$ \frac{\alpha}{\sigma_c^2} = \frac{d}{\sigma_c G_0 \tau} \tanh \left( \frac{\sigma_c G_0 \tau}{2d} \right) \quad (14) $$

Note that this equation (in $d$) has a unique solution $d(\alpha)$ if and only if $0 < \alpha/\sigma_c^2 < 1/2$. However, in view of (13), it is easier to express everything in terms of the parameter $d$ instead of $\alpha$: then the limit $\alpha/\sigma_c^2 \to 1/2$ is equivalent to $d \to +\infty$ and $\alpha \to 0$ is equivalent to $d \to 0$.

Now using the formula for computing $\sigma_Y$ we obtain by integration,

$$ \frac{\sigma_Y}{\sigma_c} = \frac{1}{2} \coth \left( \frac{\sigma_c G_0 \tau}{2d} \right) - \frac{d}{\sigma_c G_0 \tau} \quad (15) $$

To obtain the prefactor $A$ we must compute $Q^1$ which necessitates the computation of $R_{\pm}^1$. Using the continuity of $\partial_{\sigma} \mathcal{P}$ at $\pm \sigma_c$, one can find that $R_{\pm}^1(z)$ are functions satisfying

$$ \left\{ \begin{aligned} -d\partial_{z}^{2}R_{\pm}^{1} + R_{\pm}^{1} &= 0 \\ \partial_{z}R_{+}^{1}(0) &= \partial_{\sigma}Q^{0}(\sigma_{c}) = \frac{\tau G_{0}}{d\sigma_{c}\left(1 - e^{-\frac{\sigma_{c}\tau G_{0}}{d}}\right)} \\ -\partial_{z}R_{-}^{1}(0) &= \partial_{\sigma}Q^{0}(-\sigma_{c}) = \frac{\tau G_{0}e^{-\frac{\sigma_{c}\tau G_{0}}{d}}}{d\sigma_{c}\left(1 - e^{-\frac{\sigma_{c}\tau G_{0}}{d}}\right)} \end{aligned} \right. \quad (16) $$

Again, it is easy to solve this system:

$$ R_{+}^{1}(z) = \frac{\tau G_{0}}{\sqrt{d}\sigma_{c}\left(1 - e^{-\frac{\sigma_{c}\tau G_{0}}{d}}\right)} e^{\frac{-z}{\sqrt{d}}} \quad (17) $$

$$ R_{-}^{1}(z) = \frac{\tau G_{0}e^{-\frac{\sigma_{c}\tau G_{0}}{d}}}{\sqrt{d}\sigma_{c}\left(1 - e^{-\frac{\sigma_{c}\tau G_{0}}{d}}\right)} e^{\frac{-z}{\sqrt{d}}} \quad (18) $$

Now we can write down the equations satisfied by $Q^1$ using the continuity of $\mathcal{P}$ at $\sigma = \pm \sigma_c$:

$$ \left\{ \begin{aligned} -d\partial_{\sigma}^{2}Q^{1} + G_{0}\tau\partial_{\sigma}Q^{1} &= \bar{d}\left(\frac{1}{\alpha}\delta_{0} + \partial_{\sigma}^{2}Q^{0}\right) \\ Q^{1}(\sigma_{c}) &= R_{+}^{1}(0) = \frac{\tau G_{0}}{\sqrt{d}\sigma_{c}\left(1 - e^{-\frac{\sigma_{c}\tau G_{0}}{d}}\right)} \\ Q^{1}(-\sigma_{c}) &= R_{-}^{1}(0) = \frac{\tau G_{0}e^{-\frac{\sigma_{c}\tau G_{0}}{d}}}{\sqrt{d}\sigma_{c}\left(1 - e^{-\frac{\sigma_{c}\tau G_{0}}{d}}\right)} \\ \int_{-\sigma_{c}}^{\sigma_{c}} Q^{1}(\sigma)d\sigma &= 0 \end{aligned} \right. \quad (19) $$

In this system $\bar{d}$ is an unknown coefficient to be simultaneously computed with $Q^1$. The computation of $Q^1$ is tedious but straightforward. The expression of $Q^1$ are quite lengthy but
can be checked out on a symbolic computation program: if $-\sigma_c < \sigma < 0$ then

$$Q^1(\sigma) = \frac{\tau G_0}{\sigma_c \sqrt{d} \left(1 - e^{-\frac{G_0 \tau \sigma_c}{d}}\right)} e^{\frac{G_0 \tau}{d} \sigma} - \frac{\tilde{d} G_0 \tau}{d^2 \sigma_c \left(1 - e^{-\frac{G_0 \tau \sigma_c}{d}}\right)} (\sigma + \sigma_c) e^{\frac{G_0 \tau}{d} \sigma} + \frac{1}{e^{\frac{G_0 \tau}{d} \sigma_c} - e^{-\frac{G_0 \tau}{d} \sigma_c}} \times \left( \frac{\sigma_c G_0 \tau}{d} \left( \frac{\tilde{d}}{d \sigma_c} - \frac{\sqrt{d}}{\sigma_c^2} \right) + \frac{\tilde{d}}{d \sigma_c} \left( 1 + e^{-\frac{G_0 \tau}{d} \sigma_c} \right) \right) \times \left( e^{\frac{G_0 \tau}{d} (\sigma + \sigma_c)} - 1 \right) \quad (20)$$

and if $0 < \sigma < \sigma_c$

$$Q^1(\sigma) = \frac{\tau G_0}{\sigma_c \sqrt{d} \left(1 - e^{-\frac{G_0 \tau \sigma_c}{d}}\right)} e^{\frac{G_0 \tau}{d} (\sigma - \sigma_c)} + \frac{\tilde{d} G_0 \tau}{d^2 \sigma_c \left(1 - e^{-\frac{G_0 \tau \sigma_c}{d}}\right)} (\sigma - \sigma_c) e^{\frac{G_0 \tau}{d} (\sigma - \sigma_c)} + \frac{1}{e^{\frac{G_0 \tau}{d} \sigma_c} - e^{-\frac{G_0 \tau}{d} \sigma_c}} \times \left( - \frac{\sigma_c G_0 \tau}{d} \left( \frac{\tilde{d}}{d \sigma_c} - \frac{\sqrt{d}}{\sigma_c^2} \right) + \frac{\tilde{d}}{d \sigma_c} \left( 1 + e^{-\frac{G_0 \tau}{d} \sigma_c} \right) \right) \times \left( 1 - e^{\frac{G_0 \tau}{d} (\sigma - \sigma_c)} \right) \quad (21)$$

The parameter $\tilde{d}$ is selected to enforce the vanishing integral condition which amounts to taking

$$\tilde{d} = - \frac{d^{3/2} e^{\frac{\tau G_0 \sigma_c}{d}} - e^{-\frac{\tau G_0 \sigma_c}{d}} + \frac{2 \tau G_0 \sigma_c}{d}}{\sigma_c e^{\frac{\tau G_0 \sigma_c}{d}} - e^{-\frac{\tau G_0 \sigma_c}{d}} - \frac{2 \tau G_0 \sigma_c}{d}} \quad (22)$$

Finally we can compute the prefactor A which is equal to $\int \sigma Q^1(\sigma)d\sigma$. All in all, we obtain (let us note $u = (\tau G_0 \sigma_c)/d$):

$$\frac{\sigma_Y}{\sigma_c} = \frac{1}{2} \coth(u) - \frac{1}{u} \quad (23)$$

$$\frac{A}{\sqrt{\sigma_c \tau G_0}} = \frac{1}{2\sqrt{u}} \frac{e^u(u \cosh(2u) + 2 \sinh(2u) + u)}{u(\sinh(2u) - u)} - \frac{1}{2\sqrt{u}} \frac{6u^2 \coth(2u)}{u(\sinh(2u) - u)} \quad (24)$$

## Acknowledgements

KM acknowledges financial support of the French Agence Nationale de la Recherche (ANR), under grant ANR-14-CE32-0005

(project FAPRES) and FP acknowledges financial support from ERC grant ADG20110209. Most of the computations were performed using the Froggy platform of the CIMENT infrastructure (https://ciment.ujf-grenoble.fr) supported by the Rhône-Alpes region (GRANT CPER07-13, CIRA http://www.ci-ra.org/) and the Equip@Meso project (reference ANR-10-EQPX-29-01). Further we would like to thank Elisabeth Agoritsas, Eric Bertin and Jean-Louis Barrat for fruitful discussions and valuable revisions of our manuscript. Also we would like to thank Jörg Rottler for assistance in the cluster analysis of plastic events.

## References

1. P. Sollich, F. Lequeux, P. Hebraud and M. E. Cates, Phys. Rev. Lett., 1997, **78**, 2020–2023.

2. P. Hébraud and F. Lequeux, Phys. Rev. Lett., 1998, **81**, 2934–2937.

3. M. L. Falk and J. S. Langer, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 57, 7192–7205.

4. D. Rodney, A. Tanguy and D. Vandembroucq, Modell. Simul. Mater. Sci. Eng., 2011, **19**, 083001.

5. A. Nicolas, K. Martens and J. L. Barrat, EPL, 2014, **107**, 6.

6. E. Agoritsas, E. Bertin, K. Martens and J.-L. Barrat, Eur. Phys. J. E: Soft Matter Biol. Phys., 2015, **38**, 71.

7. A. S. Argon, Acta Metall., 1979, **27**, 47–58.

8. J. D. Eshelby, Proc. R. Soc. London, Ser. A, 1957, **241**, 376.

9. F. Puosi, J. Rottler and J. L. Barrat, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, **89**, 042302.

10. A. Lemaître and C. Caroli, Phys. Rev. Lett., 2009, **103**, 065501.

11. K. Martens, L. Bocquet and J.-L. Barrat, Phys. Rev. Lett., 2011, **106**, 156001.

12. K. M. Salerno and M. O. Robbins, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, **88**, 062206.

13. Z. Budrikis and S. Zapperi, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, **88**, 062403.

14. J. Lin, E. Lerner, A. Rosso and M. Wyart, Proc. Natl. Acad. Sci. U. S. A., 2014, **111**, 14382–14387.

15. L. Bocquet, A. Colin and A. Ajdari, Phys. Rev. Lett., 2009, **103**, 036001.

16. K. A. Dahmen, Y. Ben-Zion and J. T. Uhl, Nat. Phys., 2011, 7, 554–557.

17. I. Cantat and O. Pitois, J. Phys.: Condens. Matter, 2005, **17**, S3455.

18. L. Bécu, S. Manneville and A. Colin, Phys. Rev. Lett., 2006, **96**, 138302.

19. J. Goyon, A. Colin, G. Overlez, A. Ajdari and L. Bocquet, Nature, 2008, **454**, 84–87.

20. J. Goyon, A. Colin and L. Bocquet, Soft Matter, 2010, **6**, 2668–2678.

21. A. Amon, V. B. Nguyen, A. Bruand, J. Crassous and E. Clement, Phys. Rev. Lett., 2012, **108**, 135502.

22. K. Kamrin and G. Koval, Phys. Rev. Lett., 2012, **108**, 178301.

23. K. W. Desmond, P. J. Young, D. D. Chen and E. R. Weeks, Soft Matter, 2013, **9**, 3424–3436.
24 E. D. Knowlton, D. J. Pine and L. Cipelletti, *Soft Matter*, 2014, **36**, 6931–6940.

25 A. Le Bouil, A. Amon, S. McNamara and J. Crassous, *Phys. Rev. Lett.*, 2014, **112**, 246001.

26 K. E. Jensen, D. A. Weitz and F. Spaepen, *Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys.*, 2014, **90**, 042305.

27 O. Gendelman, H. G. E. Hentschel, P. K. Mishra, I. Procaccia and J. Zylberg, *Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys.*, 2014, **90**, 042315.

28 C. Liu, E. Ferrero, F. Puosi, J. Barrat and K. Martens, *arXiv:1506.08161*, 2015.

29 M. Kac, G. E. Uhlenbeck and P. C. Hemmer, *J. Math. Phys.*, 1963, **4**, 216.

30 W. Klein, H. Gould, N. Gulbahce, J. B. Rundle and K. Tiampo, *Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys.*, 2007, **75**, 031114.

31 N. Xu and C. O'Hern, *Phys. Rev. Lett.*, 2014, **94**, 055701.

32 J. Puckett and K. Daniels, *Phys. Rev. Lett.*, 2013, **110**, 058001.

33 J. Langer, *arXiv:1501.07228*, 2015.

34 S. Karmakar, E. Lerner, I. Procaccia and J. Zylberg, *Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys.*, 2010, **82**, 031301.

35 V. Mansard, A. Colin, P. Chauduri and L. Bocquet, *Soft Matter*, 2011, **7**, 5524–5527.

36 Y. Gati, PhD thesis, Ecole des ponts ParisTech, 2004.

37 E. Cances, I. Catto and Y. Gati, *SIAM J. Math. Anal.*, 2006, **37**, 60.

38 E. Cances, I. Catto, Y. Gati and C. L. Bris, *Multiscale Model. Simul.*, 2006, **4**, 1041.

39 J. Olivier, PhD thesis, Université de Grenoble, 2011.

40 J. Olivier, *Z. Angew. Math. Phys.*, 2010, **61**, 445.

41 J. Olivier and M. Renardy, *SIAM J. Appl. Math.*, 2011, **71**, 1144.

42 J. Olivier, *Sci. China-Math.*, 2012, **55**, 435–452.

43 R. Brüning, D. A. St-Onge, S. Patterson and W. Kob, *J. Phys.: Condens. Matter*, 2009, **21**, 035117.

44 T. Schneider and E. Stoll, *Phys. Rev. B: Solid State*, 1978, **17**, 1302–1322.

45 K. Salerno, C. Maloney and M. Robbins, *Phys. Rev. Lett.*, 2012, **109**, 105703.

46 A. Nicolas, J. Rottler and J.-L. Barrat, *Eur. Phys. J. E: Soft Matter Biol. Phys.*, 2014, **37**, 37–50.

47 C. Fusco, T. Albaret and A. Tanguy, *Eur. Phys. J. E: Soft Matter Biol. Phys.*, 2014, **37**, 43.

48 J. Hoshen and R. Kopelman, *Phys. Rev. B: Solid State*, 1976, **14**, 3438–3445.

49 A. Tanguy, F. Leonforte and J.-L. Barrat, *Eur. Phys. J. E: Soft Matter Biol. Phys.*, 2006, **20**, 355–364.

50 P. Sollich and A. Barra, in preparation.

51 P. Sollich, *School on Glass Formers and Glasses*, Bangalore, India, 2010.

52 H. Mizuno, S. Mossa and J.-L. Barrat, *Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys.*, 2013, **87**, 042306.