Mean-field scenario for the athermal creep dynamics of yield-stress fluids

TL;DR

This paper presents a mean-field elasto-plastic model for athermal yield-stress fluids that captures their transient creep dynamics, including re-acceleration and fluidization phases, with results depending on initial conditions.

Contribution

It introduces a novel mean-field framework that reproduces experimental creep behaviors and elucidates the physical processes during fluidization in yield-stress materials.

Findings

Reproduces non-linear shear-rate response to shear-stress.

Identifies distinct re-acceleration and fluidization phases.

Fluidization time depends on initial conditions and static yield stress.

Abstract

We develop an elasto-plastic description for the transient dynamics prior to steady flow of athermally yielding materials. Our mean-field model not only reproduces the experimentally observed non-linear time dependence of the shear-rate response to an external shear-stress, but also allows for the determination of the different physical processes involved in the onset of the re-acceleration phase after the initial critical slowing down and a distinct well defined fluidization phase. The evidenced power-law dependence of the fluidization time on the distance of the applied to an age dependent static yield stress is not universal but strongly dependent on initial conditions.

Equations

$${\partial }_t\mathcal{P}(\sigma ,t)$$

$$D(t)=\alpha \Gamma (t)\mathrm{and}\Gamma (t)=\frac{1}{\tau }{\int }_{|\sigma |>{\sigma }_c}d\sigma \mathcal{P}(\sigma ,t)$$

$$\dot{\gamma }(t)=\frac{1}{\tau G_0}{\int }_{|\sigma |>{\sigma }_c}d\sigma \sigma \mathcal{P}(\sigma ,t)$$

$${\partial }_t{\mathcal{P}}_q(\sigma ,t)$$

$${\partial }_t{\mathcal{P}}_a(\sigma ,t)$$

$$\lambda\hat{=}e^{-\big{(}\frac{s_{d}}{\sigma^{ext}}\big{)}^{2}}\;.$$

$$\frac{d}{dt}\big{(}s_{d}^{2}\big{)}=2G_{0}\sigma^{ext}\dot{\gamma}(t)-\frac{1}{\tau}\int_{|\sigma|>\sigma_{c}}d\sigma\mathcal{P}(\sigma,t)(\sigma^{2}-2\alpha)$$

$$\dot{\gamma }(t)=\frac{1}{\tau }{\int }_{|\sigma |>{\sigma }_c}d\sigma \sigma \mathcal{P}(\sigma ,t).$$

$$\frac{d}{dt}\lambda=-\frac{2G_{0}}{\sigma^{ext}}\lambda\dot{\gamma}+\frac{1}{\tau\big{(}\sigma^{ext}\big{)}^{2}}\lambda K$$

$$\frac{d}{dt}\lambda =-f(\varphi ){\lambda }^n+g(\varphi )(1-\lambda {)}^m.$$

$$\displaystyle\partial_{\sigma}\mathcal{P}\big{|}_{i}(t)=\frac{\mathcal{P}_{i+1}(t)-\mathcal{P}_{i-1}(t)}{2d\sigma}.$$

$$\displaystyle\partial_{\sigma}^{2}\mathcal{P}\big{|}_{i}(t)=\frac{\mathcal{P}_{i+1}(t)-2\mathcal{P}_{i}(t)+\mathcal{P}_{i-1}(t)}{d\sigma^{2}}.$$

$${\mathcal{P}}_i(t+dt)$$

Full text

Yield stress fluids (YSF), such as dense emulsions or pastes, display a rich rheological behavior that has attracted considerable attention recently Bonn and Denn (2009); Bonn et al. (2017). Stationary flow is typically described by a nonlinear flow curve $\sigma(\dot{\gamma})$, where $\sigma$ is the stress and $\dot{\gamma}$ the deformation rate. However, the flow curve is far from accounting for the full complexity of these systems, that involves an interplay between external driving and internal aging, potentially leading to complex thixotropic behavior. In recent years, many experiments and molecular simulations have tried to reveal this complexity using creep experiments, in which the flow rate is measured in response to a fixed stress $\sigma$ applied at a given waiting time $t_{w}$ after sample preparation Struik et al. (1978); Siebenbürger et al. (2012); Baldewa and Joshi (2012); Leocmach et al. (2014); Ballesta and Petekidis (2016); Chaudhuri and Horbach (2013); Sentjabrskaja et al. (2015); Landrum et al. (2016). These experiments, that lead, for $\sigma$ larger than a yield stress $\sigma_{Y}$, to flow or failure, reveal an intriguing behavior, with two salient features: (i) the strain-rate $\dot{\gamma}(t)$ in response to a stress larger than the yield stress is strongly non-linear and nonmonotonous, with a so called ”s-shaped” dependence of $\dot{\gamma}(t)$ Divoux et al. (2011); Siebenbürger et al. (2012); Coussot et al. (2006), including a nontrivial ”primary creep regime” often described by a power law $t^{-\mu}$. (ii) The fluidization time scale $\tau_{f}$ diverges when approaching yield stress, however in a non-universal manner.

In this work, we develop an approach that explains these features in athermal systems, in which thermal fluctuations have essentially no influence on the flow. These systems are a large subset of YSF Bonn et al. (2017), including e.g. foams, emulsions, physical gels, or granular media. In spite of the irrelevance of thermal fluctuations, the creep dynamics will depend on the initial condition determined by the preparation process and the subsequent waiting period, during which slow processes such as coarsening or compaction can alter the level of relaxation. Our approach contrasts previous attempts addressing systems in which thermal fluctuations are important, based on the soft glassy rheology model Sollich (1998); Fielding et al. (2000); Moorcroft and Fielding (2013); Merabia and Detcheverry (2016) or mode-coupling theory Fuchs and Cates (2002); Brader et al. (2009),

Our description is based on a mean-field version of the elasto-plastic scenario, that describes the flow as resulting from interactions among local plastic events triggered by the external driving, and accounts for the flow properties of athermal YSF Argon (1979); Picard et al. (2005); Nicolas and Barrat (2013); Puosi et al. (2014); Liu et al. (2016); Albaret et al. (2016); Boioli et al. (2017). It extends a previous formulation for imposed shear-rate Hébraud and Lequeux (1998) to a system subjected to an imposed stress, allowing us to address complex protocols, in particular typical creep experiments. Using this approach, we reproduce features (i) and (ii) above. We also find that the fluidization time is sensitive to initial aging, represented by an initial distribution of quenched local stresses. The divergence of the fluidization time is described by a power-law relation if the distance to yield is characterized using an age-dependent static yield stress $\sigma^{S}_{y}$. However the exponent appears, as in experiments, to be non-universal.

To investigate the creep dynamics we develop an extension of the Hébraud-Lequeux model Hébraud and Lequeux (1998) which belongs to the class of athermal, local yield stress models Hébraud and Lequeux (1998); Bocquet et al. (2009); Agoritsas et al. (2015); Agoritsas and Martens (2017). The distribution $\mathcal{P}(\sigma,t)$ characterize the stress values $\sigma$ in a subvolume of mesoscopic size. The macroscopic stress, which in creep is the externally applied load, $\sigma^{\mathrm{ext}}$ is computed as $\sigma^{\mathrm{ext}}\equiv\langle\sigma\rangle(t)\hat{=}\int\mathcal{P}(\sigma,t)\sigma d\sigma$.

The evolution equation for $\mathcal{P}$ is formulated in terms of a time dependent strain-rate, $\dot{\gamma}(t)$, as:

[TABLE]

The first term on the right hand side accounts for the local elastic response with shear modulus $G_{0}$, the second term describes local yielding at rate $1/\tau$ if the local stress exceeds a threshold $\sigma_{c}$ ($\theta$ is the Heaviside distribution) and the third term is the gain term accounting for a complete relaxation of the stress, where $\delta(\sigma)$ is the Dirac distribution. The last term is a mean-field description of the interaction generated by stress redistribution after local yielding. It describes the resulting fluctuations as a diffusive process with a time-dependent diffusion constant $D(t)$ proportional to the rate of plasticity $\Gamma(t)$:

[TABLE]

The mechanical coupling strength $\alpha$ characterizes how strongly local stresses are altered by surrounding plastic events Puosi et al. (2015). For a fixed shear-rate and $\alpha<1/2$ the model leads to a yield stress rheology, admitting a finite value of the average stress $\sigma_{y}(\alpha)<\sigma_{c}$ (dynamical yield stress) for vanishing strain-rate.

To implement a controlled stress protocol, the evolution of the shear-rate is constrained to follow the plastic shear-rate according to

[TABLE]

Using (1), it is easily verified that the average stress remains constant during the evolution of $\mathcal{P}(\sigma,t)$. The desired value $\sigma^{\mathrm{ext}}$ is imposed through the initial stress distribution $\mathcal{P}(\sigma,t=0)$.

The initial condition is defined by the stress distribution $\mathcal{P}_{0}(\sigma)=\mathcal{P}(\sigma,t=0)$. To mimic situations of a quenched system before applying the stress step, we consider distributions with zero mean $\mathcal{P}^{I}_{0}(\sigma^{\mathrm{int}})$, instantaneously shifted by the desired value of the applied stress $\sigma^{\mathrm{ext}}$ at the onset of the experiment, i.e $\mathcal{P}_{0}(\sigma)=\mathcal{P}^{I}_{0}(\sigma-\sigma^{\mathrm{ext}})$. In principle we should consider only distributions with a compact support $\sigma<\sigma_{c}$, so that the system does not evolve until the external load is applied. Hence, the model does not display aging.

In a first approach, we assume for $\mathcal{P}^{I}_{0}(\sigma^{\mathrm{int}})$ a Gaussian shape 111Strictly speaking this initial condition violates the condition of a compact support, but in practice the standard deviations studied are small enough such that the statistical weight beyond $\sigma_{c}$ can be regarded negligible. centered at zero Srolovitz et al. (1981). The only parameter is the standard deviation $s_{d}$, characterizing the level of residual heterogeneity in an amorphous system. As more relaxed systems display a less prominent Boson peak, which is indicative of a better homogeneity of the elastic properties Duval et al. (2006), we assume relaxation is also reducing the width of the stress distribution. Thus a $\mathcal{P}^{I}_{0}(\sigma^{\mathrm{int}})$ with a smaller $s_{d}$ corresponds to a more relaxed system, and we will take $s_{d}$ as an indirect measure of the age. Interestingly the standard deviation of our distribution can be formally linked to the aging parameter in the lambda-model for thixotropic materials Wei et al. (2016); Armstrong et al. (2016), as discussed in the supplementary material.

We numerically solve Eq.(1) using an explicit Euler method. The dynamic yield stress $\sigma_{y}$ of the mean-field model is a decreasing function of the mechanical coupling strength $\alpha$ Agoritsas et al. (2015); Puosi et al. (2015). Our range of investigation is restricted to $s_{d}\geq 0.22$ and $\sigma^{ext}-\sigma_{y}(\alpha)\geq 0.12$ due to numerical limitations.

*Results – * Firstly, a typical non-linear response of $\dot{\gamma}(t)$ and $\gamma(t)$ after application of a step stress $\sigma^{ext}$ larger than $\sigma_{y}$ is shown in Fig. 1. At $t\approx 1$, $\dot{\gamma}(t)$ begins to evolve significantly. For small enough $\sigma^{ext}$, the mean-field model reproduces the characteristic s-shaped curve for $\dot{\gamma}(t)$ that has been observed in various experiments Coussot et al. (2006); Siebenbürger et al. (2012); Divoux et al. (2011). For a fixed initial aging level ($s_{d}$ fixed in Fig. 1(a)), before entering the plateau of the steady-state regime, $\dot{\gamma}(t)$ displays a creep regime where the shear-rate decreases with time in an apparent power-law until it reaches a minimum. Within this creep regime, the accumulated strain $\gamma(t)$ shows a sub-linear increase in time. After the minimum in $\dot{\gamma}$, the system enters a fluidization regime where the shear-rate speeds up toward the steady-state, and correspondingly the accumulated strain $\gamma(t)$ increases super-linearly to reach the linear regime of a fluid. As the applied stress increases, the extent of the creep regime decreases, until it eventually disappears and the system enters directly the fluidization regime leading to steady-state flow (Fig. 1(a)). A similar effect of the initial aging level is found for a given applied stress( Fig. 1(b)). The duration of the creep regime decreases when increasing $s_{d}$ (decreasing age), until for a large enough $s_{d}$ the creep regime disappears and finally reaches the same stationary shear-rate. This type of dependence on the applied stress and initial aging is similar to what has been reported for bentonite suspensions Coussot et al. (2006) and colloidal hard spheres Siebenbürger et al. (2012).

Several works in the literature describe the slowing down as a power-law $\dot{\gamma}\sim t^{-\mu}$, expecting that the creep exponent $\mu$ may have universal features. Our results can indeed be fitted using such a power law, in particular when the applied stress $\sigma^{ext}$ gets small. However, we observe that the exponent $\mu$ of the apparent power-law decreases with increasing applied stress and decreasing initial aging (increasing $s_{d}$). We measure, for those $\dot{\gamma}(t)$ exhibiting a creep regime (defined as from $t=1$ to the minimum of $\dot{\gamma}$) that lasts at least one decade, the exponent $\mu$ varying from $\approx 1.2$ to $\approx 0.6$, a range comparable to the one reported in experimental studies Divoux et al. (2011); Ballesta and Petekidis (2016); Landrum et al. (2016).

We now discuss the relation between the fluidization time scale and the distance of the applied stress to the dynamical yield stress $\Delta\sigma=\sigma^{\mathrm{ext}}-\sigma_{y}$. Two time scales can be identified from the time dependence of the shear rate: $\tau_{m}$, corresponding to the minimum of $\dot{\gamma}$ and $\tau_{f}$ defined as the inflection point of $\dot{\gamma}(t)$ before entering the steady flow region. Following Divoux et al. (2011), we identify $\tau_{f}$ as the fluidization time, which can be defined even in the absence of a well developed creep regime. We will show that $\tau_{f}$ and $\tau_{m}$ are associated with different mechanisms, and $\tau_{f}$ is actually the adequate choice to characterize the time scale for entering the stationary flow.

Typical relations between $\tau_{f}$ and $\Delta\sigma$ for different aging levels are shown in Fig. 2(a). Previous experimental results Divoux et al. (2011) suggest $\tau_{f}\sim\Delta\sigma^{-\beta}$ with $\beta$ measured from $2$ to $8$ depending on the sample preparation. Note that this behavior must be distinguished from studies on thermal systems Merabia and Detcheverry (2016); Lindström et al. (2012); Sprakel et al. (2011); Gibaud et al. (2009) which suggest an exponential relation. Here our results show convexity for $s_{d}$ small (i.e. well relaxed systems) (Fig2(a)) indicating the fluidization time increases faster than a power-law as $\Delta\sigma$ approaches zero. For larger $s_{d}$, $\tau_{f}(\Delta\sigma)$ becomes closer to a power-law. The longer the waiting time before the startup of the creep protocol, the stronger the fluidization time increases for decreasing $\Delta\sigma$ (Fig. 2(a)).

To quantify the dependence of $\tau_{f}(\Delta\sigma)$ on $s_{d}$, we fit each curve in Fig 2(a) with a power law $\tau_{f}=A\Delta\sigma^{-\beta}$. As discussed above this fitting form is appropriate only for large enough $s_{d}$. For small $s_{d}$, the reported value is an effective exponent $\beta=\max\frac{d\ln\tau_{f}}{d\Delta\sigma}$, and carries qualitative information on how fast the fluidization process slows down with decreasing stress. The result is shown in the inset of Fig2(a). A systematic decrease of $\beta$ with $s_{d}$ is observed for all $\alpha$. Note that the $\beta$ value measured for carbopol microgel Divoux et al. (2011); Grenard et al. (2014) lies in the same range (from $2$ to $8$) as found here. We also notice that the exponent $\beta$ depends on the strength $\alpha$ of the mechanical coupling. A larger $\alpha$ yields a larger $\beta$ for the same $s_{d}$.

An analytic study of the complex transient dynamics is difficult, due to the strong non-linearity of Eq. (1). Therefore, we limit our analysis to a semi-qualitative discussion. Initially, the population of marginal stable nodes which is the source of plastic activity, decreases exponentially in time (the term \frac{1}{\tau}\Theta(\big{|}\sigma\big{|}-\sigma_{c})\mathcal{P}) but is also enhanced by the population of nodes having a stress slightly smaller than $\sigma_{c}$. This enhancement results from two comparable sources: a drift and a diffusion corresponding respectively to the first and second partial derivative term in Eq.(1). Note that both contributions are proportional to the population of marginally stable nodes. Further the drift and diffusion induced fluxes are, respectively, proportional to \mathcal{P}\big{|}_{\sigma=\sigma_{c}} and \partial_{\sigma}\mathcal{P}\big{|}_{\sigma=\sigma_{c}}.

Let us first consider two extreme cases with $\sigma^{\mathrm{ext}}>\sigma_{y}$, where a steadily flowing state exists according to the flow curve. If the standard deviation of the initial distribution $\mathcal{P}_{0}(\sigma)$ is large enough, the supply of unstable sites compensates the losses, and the flow will be strongly accelerated until steady state. This corresponds to the transient curves without any creep regime in Fig. 1. In contrast, if the initial Gaussian distribution has a small standard deviation, not only a small portion of the population is marginally stable, but also the values of \mathcal{P}\big{|}_{\sigma=\sigma_{c}} and \partial_{\sigma}\mathcal{P}\big{|}_{\sigma=\sigma_{c}} are close to zero. As a consequence the drift and the diffusion terms are weak and become even weaker since the marginally stable population decreases exponentially. The shear rate decreases rapidly, and the system gets stuck in a configuration where all stresses are below $\sigma_{c}$ and the flow stops. Note that this situation can be observed in experiments and simulations Bonn et al. (2017) even if the applied external stress is larger than the dynamic yield stress $\sigma_{y}(\alpha)$.

The above analysis raises the issue of the evolution of $\mathcal{P}(\sigma,t)$ that causes the transition from the creep regime to the fluidization regime and eventually the steady flow for intermediate values of $s_{d}$. It suggests further that $\tau_{f}$ diverges before $\Delta\sigma$ tends to zero (see Fig. 2(a)). For a well relaxed system (small $s_{d}$), the static yield stress $\sigma^{S}_{y}$, which is the minimum stress needed to fluidize a system at rest, is larger than the dynamical one $\sigma_{y}$. In an extreme situation where $s_{d}=0$, one should apply $\sigma^{\mathrm{ext}}\geq\sigma^{S}_{y}=\sigma_{c}$ to make the system flow. This is consistent with previous studies on transient dynamics that report on an age dependent overshoot in the stress-strain curve Fielding et al. (2000); Rottler and Robbins (2005); Rodney and Schrøder (2011). A comparison between $\sigma^{S}_{y}$ and the stress overshoot is presented in our supplementary material.

To gain a better understanding of the initial evolution of the shear-rate, we now approximate the full dynamics in the early regime, where mesoscopic subvolumes have been activated at most once, by setting $\mathcal{P}=\mathcal{P}_{q}+\mathcal{P}_{a}$ with $\mathcal{P}_{q}$ referring to the sites that have never been activated, and $\mathcal{P}_{a}$ to those activated once. Thus $\mathcal{P}_{q}(\sigma,t=0)=\mathcal{P}_{0}(\sigma)$, $\mathcal{P}_{a}(\sigma,t=0)=0$ and the distributions obey:

[TABLE]

[TABLE]

where $\dot{\gamma}_{q}(t)$, $\Gamma_{q}(t)$ and $D_{q}$ are defined as above, with $\mathcal{P}$ replaced by $\mathcal{P}_{q}$. We note that (4) and (5) approximate the full dynamics, ignoring the possibility of multiple activation. As a result they will always lead to a vanishing strain rate at long times, $\dot{\gamma}_{q}(t)\rightarrow 0$. However, the comparison between this approximation and the full solution will give us insights into the time range over which the initial condition influences the fluidization process.

The approximate solution $\dot{\gamma}_{q}(t)$ (obtained by solving (4) and (5)) and the full solution $\dot{\gamma}(t)$, for the same initial setting $(\sigma^{\mathrm{ext}},s_{d})$, are compared in Fig. 2(b). Up to the mid-fluidization, $\dot{\gamma}$ and $\dot{\gamma}_{q}$ are in good agreement, indicating that the transient dynamics is dominated by sites that undergo their first activation. Within this first stage, the flux of sites transfered from $\mathcal{P}_{q}$ to $\mathcal{P}_{a}$ first decreases with time as the population of initially unstable sites gets depleted (corresponding to the creep regime). At $\tau_{m}$, this flux starts increasing, and $\mathcal{P}_{a}$ starts to represent a significant fraction of the sites. Beyond $\tau_{f}$, $\mathcal{P}_{a}$ becomes dominant and further activation events become essential. Eventually the memory of the initial condition is lost and the steady flow regime can be achieved in the full dynamics of $\mathcal{P}$.

This analysis suggests that the critical yield stress at which fluidization takes place is not determined by the steady-state flow curve, but depends strongly on initial conditions and relaxation level. Hence, a power-law behavior cannot be expected if the dynamic yield stress is taken as reference. Rather, we estimate the static yield stress from the divergence of the fluidization time, by identifying for a given $s_{d}$ the value for which a power-law $\tau_{f}\sim\Delta\sigma_{s}^{-\beta_{s}}\equiv[\sigma^{\mathrm{ext}}-\sigma^{S}_{y}(s_{d})]^{-\beta_{s}}$ holds. Finding the best power-law fitting we estimate both $\beta_{s}(s_{d})$ and $\sigma^{S}_{y}(s_{d})$. The result of this analysis is shown in Fig. 2(c), where a power-law extends to at least one decade.

The static yield stress and exponent $\beta_{s}$ as a function of the initial aging for different couplings $\alpha$ are shown in Fig. 2(d). We use $s_{d}/(\sigma_{c}-\sigma_{y})$ as the control variable, since the effects of $s_{d}$ for different values of $\alpha$ are comparable only when measured relative to the distance between the local threshold $\sigma_{c}$ and the dynamical yield stress $\sigma_{y}(\alpha)$. When $s_{d}$ becomes of the order of $\sigma_{c}-\sigma_{y}(\alpha)$, the initial configuration contains a sufficient number of marginally stable nodes, and no difference is expected between static and dynamic yield stress. This is confirmed by the data shown in Fig. 2(d) where $\sigma^{S}_{y}-\sigma_{y}$ decreases to zero as $s_{d}/(\sigma_{c}-\sigma_{y})$ become of order $1$. In the opposite limit $s_{d}=0$, the static yield stress should be equal to $\sigma_{c}$, as discussed above. Interestingly Fig. 2(d) shows that $\sigma^{S}_{y}-\sigma_{y}$ as a function of $s_{d}/(\sigma_{c}-\sigma_{y})$ is described by a master curve independent of $\alpha$, suggesting a universal normalized relaxation level, at which the static and dynamical yield stress become identical for different systems (characterized by different $\alpha$). The insets of Fig. 2(d) shows the exponent $\beta_{s}$ against $s_{d}/(\sigma_{c}-\sigma_{y})$. A clear tendency of increasing in $\beta_{s}$ with $s_{d}/(\sigma_{c}-\sigma_{y})$ can be observed until the point where $\sigma^{S}_{y}$ becomes identical to $\sigma_{y}$. In spite of uncertainties associated with the fitting procedure, the collapse of $\beta_{s}$ for different $\alpha$ also suggests a master relation between $\beta_{s}$ and $s_{d}/(\sigma_{c}-\sigma_{y}(\alpha))$. The value of $\beta_{s}$ at large $s_{d}$ is comparable with experimental measurements Divoux et al. (2011); Grenard et al. (2014).

To summarize, we proposed a mean-field approach based on the Hébraud-Lequeux model that allows one to study various rheological protocols in athermal YSF. Using this approach , we study the creep behavior (flow induced by an imposed stress) of YSF starting from different initial levels of relaxation. We reproduce the s-shaped response of the shear-rate after imposing a step stress, in qualitative agreement with experiments Coussot et al. (2006); Siebenbürger et al. (2012). To be able to compare with experimental results, we quantify the slowing down in the creep regime with a power-law $\dot{\gamma}\sim t^{-\mu}$ and find that the exponent $\mu$ lies in the same range as reported in experiments Divoux et al. (2011). The apparent exponent is non-universal, and depends on the applied stress and on the initial relaxation level. We distinguished the different underlying mechanisms for the two time scales $\tau_{m}$ (minimum strain-rate) and $\tau_{f}$ (fluidization). $\tau_{m}$ is determined by the first plastic activations and depends sensitively on the initial distribution of internal stresses $\mathcal{P}_{0}(\sigma)$. $\tau_{f}$ characterizes the loss of memory with respect to this initial distribution. We rationalized our results for the fluidization time $\tau_{f}$ by introducing a static yield stress $\sigma^{S}_{y}$ that increases with initial relaxation. The behavior of $\tau_{f}$ can then be described by a power law $\tau_{f}\sim(\sigma^{\mathrm{ext}}-\sigma^{S}_{y})^{-\beta_{s}}$, where $\beta_{s}$ increases with decreasing initial aging, and has again values comparable to those reported in experiments. Finally, we propose that different systems can be compared by introducing an appropriate measure of relaxation based on the distribution of internal stresses, a proposal that could be tested in microscopic simulations and experiments giving access to local stresses.

*Acknowledgements – * J.-L.B. and C.L. acknowledge financial support from ERC grant ADG20110209.J-L. B is a member of Institut Universitaire de France. K.M. acknowledges financial support of the French Agence Nationale de la Recherche (ANR), under grant ANR-14-CE32-0005 (FAPRES) and CEFIPRA grant. We thank Laura Foini for detailed discussions and for her help in developing the numerical code.

Supplemental Material for:

“Mean-field scenario for the athermal creep dynamics of yield-stress fluids”

In this supplementary material we address three complementary aspects for the main article: First we explain how the static yield stress $\sigma^{S}_{y}(\alpha,s_{d})$ introduced in the manuscript is related with the stress overshoot in the stress-strain curve obtained by shear start up test in the zero shear rate limit. Second we show that an effective ”$\lambda$-model” formulation can be derived from our model. Finally we describe in detail the numerical implementation of our model used and some related numerical limits.

I Static yield stress and the stress overshoot

Stress overshoots in shear start up experiments are commonly observed in both dense hard sphere systems Rottler and Robbins (2005), but also in network forming colloidal gels Jamali et al. (2017). In the latter case, the stress overshoot is closely related with the change of structural quantities, while in the former case, which is of our interests, the stress overshoot depends on the relaxation level (age) of the system Fielding et al. (2000); Rottler and Robbins (2005); Rodney and Schrøder (2011).

The static yield stress $\sigma^{S}_{y}(s_{d},\alpha)$ is, according to our study, closely related to the initial relaxation level characterized by $s_{d}$. In this section we compare $\sigma^{S}_{y}$ with the stress overshoot during a shear start up produced by our model.

To assess the effect of initial relaxation on the stress overshoot, we perform shear start-up simulations (at constant shear rate) from an initial condition $\mathcal{P}_{0}(\sigma)$ corresponding to a centered Gaussian with different $s_{d}$, as used in the main article. The typical stress-strain curves for different $s_{d}$ at fixed $\dot{\gamma}$ and for different $\dot{\gamma}$ at a fixed $s_{d}$ are shown, respectively, in Fig. 3(a) and Fig. 3(b). The overshoot stress $\sigma_{o}$ is defined as the maximum stress of the stress-strain curve during the shear start-up. These two figures illustrate that the overshoot stress reflects the combined effects of rejuvenation by shearing and initial aging. Since the static yield stress $\sigma^{S}_{y}$ is the threshold of stress that should be applied for fluidizing a system at rest with certain initial relaxation level, one should expect that, when the shear rate becomes so slow that the rejuvenation effect becomes negligible, the overshoot stress reflecting only the initial aging should be comparable with the static yield stress.

The dependence of the overshoot stress $\sigma_{o}$ on the shear rate $\dot{\gamma}$ for a given $(\alpha,s_{d})$ can actually be fitted with a Herschel-Bulkley type relation, i.e. $\sigma_{o}=\sigma^{S}_{o}+B\dot{\gamma}^{n}$. By adjusting the value of $\sigma^{S}_{o}$, a clear linear relation can be found when $\sigma_{o}-\sigma^{S}_{o}$ is plotted against $\dot{\gamma}$ in log-log scale, Fig. 3(c). This allows us to estimate the overshoot stress in the limit of zero shear rate $\sigma^{S}_{o}$ for the given $(\alpha,s_{d})$.

In figure 3(d) we plot on top of figure(2.d) of the main article, the resulting values of $\sigma^{S}_{o}-\sigma_{y}$ against $s_{d}/(\sigma_{c}-\sigma_{y})$ for the same different values of $\alpha$. For small values of $s_{d}$, $\sigma^{S}_{y}$ agrees well with $\sigma^{S}_{o}$, which justifies our idea that the underlying physics of the static yield stress observed in creep experiments and that of the zero shear rate limit stress overshoot in shear start-up experiments are the same.

We also notice from the figure that for larger values of $s_{d}$, a small but systematic deviation exists between $\sigma^{S}_{y}$ and $\sigma^{S}_{o}$. Such initial conditions correspond to poorly relaxed systems, with a very short fluidisation time, for which deviations from the general scenario outlined in this work may be expected. In any case, a Gaussian distribution of initial stresses is unlikely to be a realistic representation of poorly relaxed systems (e.g. strongly pre-sheared systems that may rather be expected to have a non symmetric distribution of initial stresses).

II From mesoscopic mean-field description to lambda-model

We can actually derive from our model a dynamics of a specifically defined $\lambda$ fitting the general form of the lambda model Wei et al. (2016); Armstrong et al. (2016). Just as $\lambda$ does in the lambda model, $s_{d}(t)$ (the standard deviation of the distribution at time $t$) can be seen as a phenomenological parameter characterizing the aging level in our model. However, as the model discusses the evolution of the full distribution $\mathcal{P}(\sigma,t)$, it is richer than a model that only uses one scalar parameter. Still, it is possible to make an approximate mapping of $s_{d}(t)$ to $\lambda(t)$ by defining:

[TABLE]

As so, when $s_{d}$ gets small, $\lambda$ gets close to one corresponding to fully aged system, and when $s_{d}$ gets large, $\lambda$ gets close to zero corresponding to a rejuvenated system Wei et al. (2016); Armstrong et al. (2016).

The time evolution of $\lambda$ is then determined by that of $s_{d}^{2}=\langle\sigma^{2}\rangle-\langle\sigma\rangle^{2}=\langle\sigma^{2}\rangle-(\sigma^{ext})^{2}$. Since $\sigma^{ext}$ is fixed, \frac{d}{dt}\big{(}s_{d}^{2}\big{)}=\frac{d}{dt}\langle\sigma^{2}\rangle. By integrating on both sides of Eq(1) in the main text of paper with $\int d\sigma\sigma^{2}\bullet$, one obtains:

[TABLE]

with

[TABLE]

Here for simplicity, we only consider a positive stress is applied, which insures $\dot{\gamma}(t)\geq 0$. Since $\sigma_{c}=1$ and $\alpha<0.5$, the integration in the second term of Eq.7 is strictly positive. This also gains clear interpretation: the first term in Eq.7 representing the flow contributes to rejuvenate the system, while the second term representing the local relaxations contribute to aging. By changing variables, one easily obtains the dynamics of $\lambda$:

[TABLE]

with $K(>0)$ is the integration in the second term of r.h.s. of Eq.7. $\dot{\gamma}$ and $K$ are in effect two functionals of $\mathcal{P}(\sigma,t)$, however they can be roughly regarded as two functions of $\lambda$ and $\sigma^{ext}$ if a simple parametrization (Gaussian or exponential) of $\mathcal{P}(\sigma,t)$ is chosen.

Comparing Eq.9 with the general form of lambda model Wei et al. (2016):

[TABLE]

Both formulations contain, on the r.h.s., a positive aging term and a negative term of rejuvenation driven by the external loading. It is not difficult to verify that the aging term in Eq.9 vanishes at $\lambda\rightarrow 1$ and the rejuvenation term vanishes at $\lambda\rightarrow 0$. Both the aging and rejuvenation terms are positive functions of $\sigma^{ext}$. Besides keeping the consistency with the lambda model, Eq.9 not only introduces more complexity on the interplay between the external loading and the structural parameter $\lambda$, but also offers a more generic physical picture of the widely studied lambda model from a mesoscopic point of view, i.e. the structural parameter can be related with mesoscopic stress fluctuations.

III Numerical Implementation

We adopt the finite-difference Euler integration method to numerically solve the mean-field model described by Eq.(1) and the schematic evolution described by Eq.(4) and Eq.(5) in the main article.

We discretize the support of the time dependent probability density of local shear stresses $\mathcal{P}(\sigma,t)$ by $d\sigma$ within a domain of $\left[\sigma_{\mathrm{min}},\sigma_{\mathrm{max}}\right]$ and discretize the time by $dt$. The total number of points in stress space is then $N=(\sigma_{\mathrm{max}}-\sigma_{\mathrm{min}})/d\sigma$. The probability density $\mathcal{P}(\sigma,t)$ is therefore discretized into $\mathcal{P}_{i}(t)=\mathcal{P}(\sigma_{i},t)$, with $\sigma_{i}=\sigma_{\mathrm{min}}+i\times d\sigma$ and $i=0,1,...,N-1$. A periodic boundary condition is used, so that $P_{N}=P_{0}$. The source term Dirac function $\delta(\sigma)$ is approximated by a discretized Gaussian function with a narrow standard deviation, i.e. $\delta_{i}=\frac{1}{\xi_{1}\sqrt[]{2\pi}}\exp(-\sigma_{i}^{2}/\xi_{1}^{2})$ and the step function $\theta(|\sigma|-\sigma_{c})$ is approximated by $\theta_{i}=\left[\tanh(-\xi_{2}(\sigma_{i}+\sigma_{c}))+\tanh(\xi_{2}(\sigma_{i}-\sigma_{c}))+2\right]/2$. The two parameters $\xi_{1}$ and $\xi_{2}$ control respectively the width of the Gaussian and the width of the region with large variation of $\theta_{i}$. Taking the limits $\xi_{1}\rightarrow 0$ and $\xi_{2}\rightarrow\infty$, $\delta_{i}$ and $\theta_{i}$ converge to $\delta(\sigma)$ and $\theta(|\sigma|-\sigma_{c})$ respectively, along which we should also take $d\sigma\rightarrow 0$ for the numerical discretization reflects well the regularity of these functions, whereas the trade off being the computing cost.

We use a symmetrised finite difference method for estimating the first and second partial derivatives with respect to stress, as follows:

[TABLE]

[TABLE]

Periodic boundary conditions are applied for estimating the partial derivatives at the two edges. We therefore compute at each time step $t$ with the current probability density $\mathcal{P}_{i}(t)$ all necessary quantities for performing a Euler integration for the time evolution, i.e.

[TABLE]

where $\Gamma(t)=\frac{1}{\tau}\sum_{i}d\sigma\theta_{i}\mathcal{P}_{i}(t)$. For implementing the shear rate control protocol, $\dot{\gamma}(t)$ is simply set to a constant equal to the desired value. For implementing the stress control protocol, we compute $\dot{\gamma}(t)=\frac{1}{G_{0}\tau}\sum_{i}d\sigma\theta_{i}\sigma_{i}\mathcal{P}_{i}(t)$.

We compared the solution of a pure diffusion equation solved by our numerical method and the analytical solution with initial conditions as a narrow Gaussian centered at zero. Generally the condition $dt\lesssim d\sigma^{2}/2$ is required for the numerical solutions to be stable and to approximate well the true solution for a diffusion equation. As long as the width of the probability density is much less than the numerical domain $[\sigma_{\mathrm{min}},\sigma_{\mathrm{max}}]$, our numerical solution compares well with the analytical one (data not shown here).

From the analytical stationary solution of the Hébraud-Lequeux model Agoritsas et al. (2015); EliArxiv16 , we know that, for the stress range that we are interested for studying the creep behavior, the probability density $\mathcal{P}(\sigma)$ is mainly weighted between [math] and $\sigma_{c}=1$ and the stationary solution beyond an absolute stress value of two $\mathcal{P}_{stat}(|\sigma|>2)$ is quasi-null. Thus we restrict the numerical domain constraint by $\sigma_{\mathrm{min}}=-5$ and $\sigma_{\mathrm{max}}=5$, which is justified by our numerical solutions with $\mathcal{P}_{i=0}\approx\mathcal{P}_{i=N-1}\lesssim 10^{-8}$ during its evolution. As mentioned before, the stress discretization $d\sigma$ is determined by how much we want $\theta_{i}$ and $\delta_{i}$ to be, respectively close to the real step function and the Dirac function, so that we try to take as small as possible the numerical domain to reduce the computational cost with respect to the number of points in stress space.

Our discretization $d\sigma$ is chosen according to the following considerations. Firstly, the spacial precision around $\sigma=\pm\sigma_{c}=\pm 1$ is crucial for well solving the dynamics for some parameter settings when $s_{d}$ and $\sigma^{ext}$ are both small, because a very small part of the probability density is present beyond $\pm\sigma_{c}$. It has been shown that the dynamics of this boundary layer is crucially determining the overall dynamics EliArxiv16. Another constraint is coming from the regularization of the derivatives of the Heaviside function. The width of the approximated function must be small enough with respect to the width of the boundary layer (given by $\mathcal{P}(|\sigma|>\sigma_{c})$) to approximate well the dynamics. Thus the choice of $d\sigma$ needs to be small enough to both well discretize the regular form of $\mathcal{P}(\sigma)$ and the width of the approximated step function $\theta_{i}$, otherwise some uncontrollable errors may occur. But note that decreasing $d\sigma$ increases the computing cost within one time step. Besides, since $dt\lesssim d\sigma^{2}/2$ is required, fluidization time being fixed, decreasing $d\sigma$ also means increasing the number of computing time steps for reaching the stationary state during one creep simulation. Taking all the above constraints into account, we have to tune $d\sigma$ and $dt$ to compromise between the computing cost and the domain of parameter settings $(\sigma^{ext},s_{d})$ that can be faithfully explored by our code. For obtaining the results in the main article, we set $dt=10^{-5}$ and $d\sigma=(\sigma_{\mathrm{max}}-\sigma_{\mathrm{min}})/8192\approx 10^{-3}$, which is at the limit of the stability condition $dt\lesssim d\sigma^{2}/2$. The results coming from this numerical discretization are compared with those coming from a numerical discretization where $dt$ fully satisfying the stability condition, and only a tiny difference is observed which does not affect quantitatively the physics of creep. As a benchmark of our numerical method, we compare our numerical steady state flow curves with those obtained by implicit analytical stationary solutions of the Hébraud-Lequeux model Agoritsas et al. (2015); EliArxiv16. The good agreement shown in Fig.4(a) validates our simulations.

With the above discretization settings and within the acceptable simulation duration, for some parameter settings with both small values of $s_{d}$ and $\sigma^{ext}$, we still did not observe the fluidization regime in the creep curve $\dot{\gamma}(t)$ and the shear rate $\dot{\gamma}$ goes down with time until the end of simulation, as example in Fig.4(b). With our hypothetical relation \tau_{f}=A(s_{d},\alpha)\big{[}\sigma-\sigma^{S}_{y}(s_{d},\alpha)\big{]}^{-\beta_{s}(s_{d},\alpha)} justified by the figure(2.c) in the main article, we extrapolate the $\tau_{f}(s_{d},\sigma;\alpha)$ for these parameter settings if $\sigma>\sigma^{S}_{y}(s_{d},\alpha)$. If $\sigma<\sigma^{S}_{y}(s_{d},\alpha)$, $\tau_{f}=\infty$ meaning that at some point $\dot{\gamma}(t)$ vanishes exponentially. It is found out that in the first case when $\sigma>\sigma^{S}_{y}(s_{d},\alpha)$, the $\tau_{f}$ extrapolated for these parameters are out of the scope of our total simulation time and in the second case $\sigma<\sigma^{S}_{y}(s_{d},\alpha)$, $\dot{\gamma}(t)$ decreases faster than a power law, as the bottom curves in Fig. 4(b). These observations are not contradictory with but even support our interpretations.