Creep dynamics of athermal amorphous materials: a mesoscopic approach

Yield stress fluids display complex dynamics, in particular when driven into the transient regime between the solid and the flowing state. Inspired by creep experiments on dense amorphous materials, we implement mesocale elasto-plastic descriptions to analyze such transient dynamics in athermal systems. Both our mean-field and space-dependent approaches consistently reproduce the typical experimental strain rate responses to different applied steps in stress. Moreover, they allow us to understand basic processes involved in the strain rate slowing down (creep) and the strain rate acceleration (fluidization) phases. The fluidization time increases in a power-law fashion as the applied external stress approaches a static yield stress. This stress value is related to the stress over-shoot in shear start-up experiments, and it is known to depend on sample preparation and age. By calculating correlations of the accumulated plasticity in the spatially resolved model, we reveal different modes of cooperative motion during the creep dynamics.

Yield-stress fluids (YSFs), such as dense emulsions, colloidal suspensions and pastes, display a rich rheological behavior that has attracted considerable attention in the last decades (see reviews by Bonn et al. [1] and Nicolas et al. [2]). Typically, the rheological behavior of YSFs is characterized by the flow curve measured in a stationary flowing state. The dependency of the stationary shear stress Σ on the applied shear rateγ is in many cases well described by a generalised form Σ(γ) ≈ σ Y + Aγ n , where the prefactor A, the "Herschel-Bulkley" exponent n and the yield stress σ Y (better referred as the dynamical yield stress) are the relevant fitting parameters. But it is well known that assessing the material's bulk properties at finite shear flow and in the steady state limit, does not fully account for the complex dynamics of certain YSFs. For example, the interplay between external driving and internal aging has been shown to cause complex thixotropic behaviors [3], leading to non-homogeneous flow even under homogeneous driving conditions. The important challenge is then to study not only the well established flow properties in the homogeneous steady flow regime, but also the spatially resolved transient dynamics that bridges the solid response at small deformations to the flowing state at large deformations.
In recent literature, efforts have involved typically two kinds of protocol: shear start-up and creep tests [4][5][6][7]. The first approach controls the applied strain whereas the latter records the strain rate response to a sudden stress step. In particular, creep experiments measure the strain rate evolution in response to a fixed stress σ applied at a given waiting time t w after sample preparation [4,5,[8][9][10][11][12][13]. In this way, the response of the system is probed as a function of its initial age. Such experiments reveal an intriguing behavior with two salient features: (i) the strain-rateγ(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γ(t) [5][6][7], including a nontrivial "primary creep regime" often de-scribed by a power lawγ ∼ t −µ ; (ii) the fluidization time scale τ f diverges when approaching the yield stress, yet in a non-universal manner. Both features are found to depend on sample preparation. The experiments ultimately lead for small applied stresses to a dynamical arrest or steady creep and for sufficiently large applied stresses to steady flow or failure, depending on the material.
In this work, we reproduce and interpret the above experimental features using mesoscopic modeling approaches. The implementations we use are suitable for athermal amorphous systems, which constitute a large sub-class of YSFs, including foams, emulsions, physical gels and granular media [1]; where large Peclet numbers assure that thermal fluctuations are negligible compared with mechanical fluctuations induced by the response to an external driving. Note that the material mechanical properties and dynamics will always depend on its preparation protocol and subsequent waiting period prior to deformation; during which slow temperature-dependent processes such as glassy relaxation are indeed relevant. In any case, our approach complements previous studies addressing creep in systems for which thermal fluctuations are also important during the dynamics, based on the soft glassy rheology model [14][15][16][17] or mode-coupling theory [18,19]. We work in particular with mean-field [20] and spatially-resolved versions of elasto-plastic models; where the flow of YSFs results from the interaction among local plastic events triggered by external driving. It has been shown that such models account for various properties of YSFs [21][22][23][24][25][26][27]. While previous formulations consider a fixed strain rate as a control parameter [25,[28][29][30][31], we extend those models here to describe the constant imposed stress condition relevant for creep. With that, we reproduce the experimental features (i) and (ii) explained above. Interestingly, we find the same fluidization time dependence on initial aging as Siebenbürger et al. [5]. Furthermore, we show that the divergence of the fluidization time is described by a power-law relation when the distance to yield is measured with respect to an agedependent static yield stress σ s Y larger than the dynamical yield stress σ Y itself. As in experiments, the resulting power-law exponent appears to be non-universal.

Mesoscopic elastoplastic model
Our spatially-resolved approach is extended from a previous version used to describe steady state flows [25,28,29,32]. It consists in a square lattice where a node (ij) represents a typical cluster of particles undergoing plastic events [21] and i, j the discretized coordinates along the x and y directions respectively. Once a reference state is set, each node is associated with a local plastic deformation γ pl ij which is, in general, heterogeneous. In addition, a node can develop an elastic strain γ el ij associated with a local stress σ ij = µγ el ij . The local stress is composed of two parts, where σ EXT is the externally applied uniform stress, and σ IN T ij encodes the stress fluctuations resulting from the interactions between plastic heterogeneities, more precisely The interaction kernel G is of the Eshelby's type [33] as described in [29], plus an homogeneous part 1/N with N the system size, so that the integral over space of the internal stress caused by any plastic strain field is null. Thus σ IN T ij describes the local stress fluctuations in a macroscopically stress free state. Applying a macroscopic stress amounts to shifting uniformly the local stress without altering internal fluctuations.
Besides, each node alternates between a local plastic state and a local elastic state by switching a local state variable n ij between 1 and 0,respectively. Explicitly when |σ ij | is larger than σ c , the site becomes plastic (n ij = 0 −→ n ij = 1) at a rate 1/τ pl and becomes elastic again (n ij = 1 −→ n ij = 0) at a rate 1/τ res [22,29], i.e.
The local dynamics is expressed as We set in our simulationτ = τ res = τ pl = 1, σ c = 1 and µ = 1. The model described above is a reformulation of the model in [29] but allowing us to model a stress controlled protocol. To summarize, Eq. (1) and Eq. (3) forms a closed stochastic dynamical system governing the evolution of the plastic strain field γ pl ij (t). Given the initial condition defined by γ pl ij (t = 0) and the imposed stress σ EXT , we simulate the system and measure the macroscopic plastic strain γ pl (t)= ij γ pl ij (t)/N , from which the strain rate response γ pl (t) results directly.

Mean-field approach for creep dynamics
The probability distribution function P (σ, t)dσ gives the fraction of nodes with a local stress in [σ, σ + dσ] at time t. Our mean-field approach to approximating the time evolution of this probability distribution for a typical site is inspired by the Hébraud-Lequeux model [30] and thus belongs to the class of athermal local yield stress models [31,34]. One should note that for the derivation of this model, we assume several strong simplifications with respect to the spatial elasto-plastic description, such that all comparisons can only be done on a qualitative level. As we will discuss later the detailed aspects of the creep curve depend on these simplifications whereas other more general features, like power-law scalings, appear to be very general. Our mean-field approach describes the dynamics of the distribution P (σ, t), and differs from the original model [30] by further taking into account a strain rate that may vary in timeγ(t): The first term on the right hand side describes local yielding with rate 1/τ if the local stress exceeds a yield stress σ c (θ denotes the usual Heaviside distribution) and the second term is the corresponding gain term accounting for a instantaneous complete relaxation of the stress.
Here δ(σ) is the Dirac distribution, and Γ(t) is the rate of plasticity The third term accounts for the local elastic response with shear modulus G 0 . The last term in Eq. (4) is a mean-field approximation for interactions between different macroscopic regions presented in the spatially resolved model. This term describes in an effective way the stress fluctuations caused by the elastic response to surrounding yielding events as a diffusive stochastic process with a time-dependent diffusion constant D(t) proportional to the rate of plasticity Γ(t): This approximation by a diffusive process assumes that the "kicks" received by a typical site are uncorrelated and of finite variance, an approximation to the true dynamics of elastoplastic models that has been improved, in the quasistatic limit of vanishing strain rate, by the work of Lin and Wyart [35] taking into account the fact that these kicks have a broad distribution. We are not aware of any similar extension to finite strain rates, and in this work will remain at the simpler level of description.
Driving the system with a constant stress, dσσP (σ, t) = constant, leads by integration of 4 to a self consistent determination of the strain rate aṡ Note that this strain rate corresponds physically to the rate of total plastic strain released by the yielding sites, each site releasing σ/G 0 , Equations 4 and 7 can be solved numerically starting from a given initial condition P 0 = P (σ, t = 0).

Initial condition and aging
In the mean-field approach, the initial condition is fully defined by P 0 (σ) = P (σ, t = 0). To represent the quenched state of a system before applying the step stress, we consider a distribution of zero mean P I 0 (σ) that describes the internal stress fluctuations of a system in a macroscopically stress free state. This distribution will be instantaneously shifted to the desired value of the imposed stress σ EXT at the beginning of the creep protocol to mimic the application of a stress step, i.e. P 0 = P I 0 (σ − σ EXT ). In the spatially resolved elastoplastic model, the initial state is defined by γ pl ij (t = 0). Using Eq. (1), γ pl ij (t = 0) can be converted to a field of internal stress fluctuations σ IN T ij (t = 0) from which a zero mean distribution, such as P I 0 (σ IN T ), can be constructed. Then the distribution of local stresses at the onset of creep experiments can also be described as an instantaneous shift of the zero mean distribution with the imposed step stress σ EXT , such as P I 0 (σ − σ EXT ). In practice, we first choose a specific form of P I 0 (σ), then convert it back to a random realization of γ pl ij (t = 0) and the creep test is simulated by evolving the model under a fixed value of σ EXT . Once an initial condition is prepared, we numerically integrate the dynamic equations using an explicit Euler method.
In principle we should consider only distributions P I 0 (σ) with a compact support in both models, i.e. P I 0 (σ > σ c ) = 0, so that the system does not evolve until the external load is applied. Hence, our models do not explicitly resolve the aging dynamics, but we mimic the role of aging by using different choices of the initial condition. In a first approach, we assume for the distribution P I 0 (σ int ) a Gaussian shape 1 centered at zero [36]. The only parameter is the standard deviation s d , characterizing the level of residual heterogeneity in the amorphous system. As more relaxed systems display a less prominent "Boson peak", indicative of a better homogeneity of the elastic properties [37], we assume that relaxation is also reducing the width of the stress distribution. Thus a P I 0 (σ int ) with a smaller s d corresponds to a more relaxed system, and we take 1/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 [38,39], as discussed in the Supplementary Material of Ref. [20].

Flow curves
Let us start, prior to the study of the creep behavior, by comparing the flow curves obtained from the stress control protocol here presented to those obtained from the strain rate control protocol explained elsewhere [29,30]. In the stress control protocol, we set the stress by choosing an arbitrary initial condition among those with the desired stress value. The system evolves with a stresspreserved dynamics and reaches a steady flow regime at large times when the memory of the initial condition is completely erased. We then averageγ(t) over a large time window in the steady state. For both the spatial model and the mean-field approach, the comparison shown in Fig. 1 reveals a good consistency between the two types of protocols. The dynamic yield stress σ Y of the spatiallyresolved model is estimated to be ∼ 0.7536. The dynamic support, but in practice the standard deviations studied are small enough such that the statistical weight beyond σc can be regarded as negligible.
yield stress σ Y of the mean-field model is a decreasing function of the mechanical coupling strength α, as explained in Ref. [31].

Creep curves
In the following, we use exclusively the stresscontrolled protocol. An initial condition with an imposed stress is chosen, corresponding to the application of a step stress on a relaxed sample in experiments. Typical responses ofγ(t) and γ(t) for the two models just after the application of a step stress are shown in Fig. 2 and Fig. 3. All these curves are obtained with an imposed stress larger than the dynamical yield stress, i.e. ∆σ ≡ σ EXT −σ Y > 0. The fact that the two models differ in their behaviors ofγ(t) for t 1 is due to the different ways in which they describe the plasticity, namely a viscous relaxation in the spatially resolved model compared to an instantaneous one in the mean field approach. Beyond a microscopic time scale t mic 1, up to which they depend on the details of the different dynamics, we notice that the strain responses of the two models behave qualitatively in the same way. Already at early times after t mic 1, the strain rateγ monotonically decreases for the smallest ∆σ > 0 cases and correspondingly the strain γ reaches a plateau. This implies that even when the applied stress step is above the dynamic yield stress, a quenched system submitted to a creep test may not yield to a flowing state. This is one of the main differences with fixed shear rate protocols, where the system is always forced to yield. We will come back to this point later. For a larger imposed stress, instead, both the spatial and mean-field models reproduce the characteristic S-shaped curve forγ(t) observed in experiments [5][6][7]. For a fixed initial aging level (s d fixed in Fig. 2), before entering the plateau of the steady-state regime, γ(t) displays a creep regime where the strain rate decreases with time in an apparent power-law fashion until it reaches a minimum. Within this creep regime, the accumulated strain γ(t) shows a sub-linear increase in time. After the minimum inγ, the system enters a fluidization regime where the strain rate speeds up toward the steady-state, and correspondingly the accumulated strain γ(t) increases super-linearly to reach the linear regime of a steady flowing state. As we further increase the stress, the extent of the creep regime decreases, until it eventually disappears and the system enters directly the fluidization regime after t mic .
A similar effect is found for a given applied stress when we vary the initial aging level (Fig. 3). The duration of the creep regime decreases when increasing s d (decreasing age), up to the point where it disappears for large enough s d , and all curves reach the same steady strain rate. This indicates that, at a given applied stress, a less relaxed system is more likely to be fluidized. These dependences on the applied stress and the initial relaxation are reminiscent of the creep tests performed in bentonite suspensions [7] and colloidal hard sphere systems [5].
Several works suggest a power law slowing downγ ∼ t −µ in the creep regime [6,11,13,40,41]. Our data can be indeed fitted with such a power law for a modest range of the parameters (s d , σ EXT ). By doing so in curves where we can fit at least one decade of power-law, we observe that the exponent µ decreases with increasing applied stress and decreasing initial aging (increasing s d ). The values of µ extracted from the two models are similar and vary from 0.6 to 1.2 (fits not shown here) depending on s d and σ EXT , a range comparable to those reported in experiments [6,11,13].
Note that, in order to produce a comparable time dependence of the strain rate, significantly different values of the control parameters ∆σ and s d are used as input in the two different models, as shown in Fig. 2 and Fig. 3. In principle some of the differences in the creep response could result from finite size effects in the spatial model and discretization effects in the mean-field description. But here, we took care that for the parameter range studied, these effects are not relevant. Instead, to understand this discrepancy, one has to recall that the rules for the local plastic deformation in the two models are very different. In the mean-field model the local stress release is instantaneous and complete whereas the spatial model implements a duration of events and only a partial local stress relaxation. This is also the origin of different behaviors of the two models in theirγ(t) for t t mic 1. Before the vertical dashed line ( Fig.2 and Fig.3), the strain rateγ(t) from the elasto-plastic model increases linearly (thus γ(t) increases quadratically) until t ≈ t mic , while the strain rate from the mean-field model begin with a finite value and varying significantly only after t mic 1. Another obvious difference is the interaction kernel of the spatial model, leading to a spatially correlated dynamics. As shown before, the flow curves do not match and in the mean-field model, they depend on the value of α, and consequently so do the creep curves. For these reasons we can only hope to qualitatively reproduce the creep curves with the mean-field model. Interestingly, however there are other quantities like the fluidization time dependence on the imposed stress that are still quantitatively comparable and thus point towards more general behavior. We will discuss this point in detail in the following section.

Fluidization time
We now discuss the relation between the fluidization time scale and the distance of the applied stress to the dynamical yield stress ∆σ = σ EXT −σ Y . Two characteristic times can be recognized in the evolution ofγ, at least for intermediate values of the applied stress. We call τ m the time where the minimum ofγ(t) occurs, after the creep regime, and define τ f as the inflection point ofγ(t) in the fluidization regime before entering the steady state. Following [6], we choose τ f to characterize a fluidization time scale. Note that τ f is always defined, even in the absence of a creep regime (e.g., for large values of ∆σ) where we can still recognize a fluidization phase with an acceleration ofγ, and it will characterize the typical time needed to reach the steady flow.   4 shows the τ f dependence on ∆σ for different initial aging levels. Experimental results on a carbopol microgel [6] suggest a power law dependence τ f ∼ ∆σ −β with β measured from 2 to 8 depending on sample preparation. This behavior must be distinguished from studies on thermal systems [17,[42][43][44] which suggest an exponential relation instead. Our results show, especially when s d is small (well relaxed systems), a convexity of the curves indicating that the fluidization time increases faster than a power law as ∆σ approaches zero. As s d becomes larger, τ f (∆σ) becomes more power-law-like. Fig. 4 also shows that more relaxed systems display a stronger increase in τ f for decreasing ∆σ. This aging dependence of τ f (∆σ) qualitatively agrees with the experimental results of Divoux et al. [6].
Despite the different functional forms when changing s d , to quantify the dependence of τ f (∆σ) on aging, we define for each curve a value of β. When appropriate (large s d , as discussed above) we directly fit a power law τ f = A∆σ −β , finding in both the spatial and mean-field models consistent values (β ≈ 2.2 and 2.38, respectively), lying in the range of experimental results. For initially well relaxed systems (small s d ), we define an effective β = max d ln τ f d∆σ , that equally carries information on how fast the fluidization process slows down with decreasing stress. Results are shown in Fig. 5(a). We see a systematic decrease of β with s d for both models. Moreover, the estimated values of β cover the same range (from 2 to 8) found in carbopol microgels [6,45]. We also notice that, for the mean field model the exponent β depends on the mechanical coupling strength α. A larger mechanical coupling strength yields a larger value of β for the same s d . It has been proposed in earlier works [46] that the value of α should depend on the specific form of the interaction kernel, and within simulations of particle based models it has been found that α typically lies in the range [47] of 0.26 to 0.33, which is consistent with the range of α of the mean-field model studied here.

Rationalization of the creep dynamics
In this section, we attempt to rationalize the behavior of the global observables reported above by analyzing in detail the evolution of the stress probability distribution. Although the two models seem quite different in their formulation, the underlying physical process are quite similar. For example, local plasticity in the mean-field model consists of a total release of local stress and a sudden return to the elastic state, which can be viewed as the local plasticity in the elasto-plastic model in the limit τ res → 0 andτ → 0. In the elasto-plastic model, the stress released by plastic events is re-distributed according to the interaction kernel through Eq. (1) keeping the average stress constant but broadening the distribution P (σ, t) with an amplitude roughly proportional to the rate of plastic events. This effect of stress redistribution is reasonably approximated in the mean-field model by the diffusive term in Eq. (4) with the diffusion coefficient proportional to the plastic activation rate. Although the more realistic elasto-plastic model contains more infor- mation, the simpler mean-field model already captures the key ingredients needed to understand the underlying mechanism of creep in athermal amorphous systems.
We therefore focus on the mean-field model and, to gain insights into the complex dynamics during creep tests, proceed to analyse the complex nonlinear behaviour of Eq. (4) in a qualitative manner. The time evolution of P (σ, t) is driven by the existence of a population of sites in an unstable state, (P (σ > σ c )). Actually, by integrating Eq. (4) beyond σ c , one obtains using that the negative part beyond −σ c is negligible for positively imposed stress andγ(t) ≈ σc G0 Γ(t), (since P (σ) weights little and decreases fast beyond σ c ). When the population of unstable sites τ Γ is non-zero, it decreases exponentially due to the loss term of plastic activation represented by the first term on the r.h.s of Eq. (8). At the same time, it is supplied by two comparable fluxes represented by the last term in Eq. (8): the stress drift due to the elasticity and the stress diffusion due to stress redistribution. These fluxes are both proportional to the unstable population. In particular, the drift and diffusion induced fluxes are respectively proportional to P σ=σc and ∂ σ P σ=σc .
Let us first consider two extreme cases with σ EXT > σ Y , where a steady flowing state exists according to the flow curve. If the standard deviation of the initial Gaussian P 0 (σ) is large enough, the supply of the unstable sites is as important as the plastic activation. Thus, the supply and the loss rapidly reach a situation where they compensate each other. This corresponds to the curves with no creep regime at high imposed stresses in Figs. 2 and 3. On the other hand, if the standard deviation of the initial distribution P 0 (σ) is very small, not only a small portion of the population is unstable but also the values of P σ=σc and ∂ σ P σ=σc are close to zero. The term of supply is then negligible compared to the loss term in Eq. (8). As a consequence, drift and diffusion are as weak as the unstable population at the beginning and become even weaker as the unstable population decreases exponentially. The strain rate rapidly decreases, the system gets eventually stuck in a configuration where all sites are below σ c and the flow stops. This situation corresponds to the curves with vanishing strain rate in Fig. 2. Note that this situation can be observed in experiments and simulations with fixed stress protocols [1] even if the applied external stress is larger than the dynamic yield stress σ Y (α).
The above analysis raises the question of the evolution of P (σ, t) that causes the transition from the creep regime to the fluidization regime and eventually the steady flow, observed for intermediate values of s d . It further suggests that τ f can diverge before ∆σ tends to zero (see Fig. 4).
For a well relaxed system (small s d ), we can therefore introduce a static yield stress σ s Y defined as the minimum stress needed to fluidize a system at rest. σ s Y will depend on the initial state and will be larger than the dynamical yield stress σ Y . In an extreme situation where s d = 0, one should apply σ EXT ≥ σ s Y = σ c to make the system flow. This is consistent with previous studies on transient dynamics that report an age dependent overshoot in the stress-strain curve [15,[48][49][50]. Actually, the stress overshoot in the zero strain rate limit is closely related with σ s Y . A comparison between the two is discussed later on. To gain a better understanding of the initial evolution of the strain rate, we study the full dynamics in the early regime -such that mesoscopic blocks have been activated at most once-by setting P = P q +P a , with P q referring to the sites that have never been activated and P a to those activated once. Thus P q (σ, t = 0) = P 0 (σ), P a (σ, t = 0) = 0 and the distributions obey: whereγ q (t), Γ q (t) and D q are defined as above, with P replaced by P q . We note that Eq. (9) and Eq. (10) 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,γ q (t) → 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.
In Fig. 6, we compare the approximate solutionγ q (t) obtained by solving Eqs.9 and 10 with the full solutioṅ γ(t), for the same initial settings (s d , σ EXT ). We used three different values of the applied stress to have different extents of the creep regime. Fig. 6 confirms that in all situations, up to the mid-fluidization regime,γ q (t) anḋ γ(t) are in good agreement, indicating that the dynamics governing the creep regime is dominated by the sites undergoing their very first activations. The curves also show that multiple plastic activations must come into play for the crossover from the fluidization regime to the steady flow to take place. We can conclude that the existence of a creep regime and the value of τ m is determined by the initial condition, while the full fluidization process, characterized by τ f , corresponds to a process of diluting the memory of the initial state through multiple plastic activations. The previous analysis suggests that the critical yield stress above which the fluidization takes place, is not determined by the flow curve, but rather by the initial state and its level of relaxation. Hence, the divergence of τ f as a power-law cannot be expected if the dynamic yield stress σ Y is taken as a reference. We then estimate a static yield stress from the divergence of the fluidization time τ f , by identifying for each s d the value σ s Y (s d ) for which a power-law τ f ∼ ∆σ −βs s ≡ [σ ext − σ s Y (s d )] −βs holds. Finding the best power-law fitting we estimate both β s (s d ) and σ s Y (s d ). The result of this analysis is shown in Fig. 7. We observe power-laws that span over a decade each, with effective exponents depending on s d .
The exponent values display a clearer trend in the meanfield case, where β s seems to increase systematically with s d . We note however that, for the elasto-plastic model, finite size effects may delay a bit the fluidization τ f when the imposed stress is small, so that the static yield stress may be overestimated 2 .
In order to make more explicit the dependency of σ s Y and β s on s d , we plot σ s Y −σ Y and β s against s d /(σ c −σ y ) in Figs. 8 and 9 (circles). This choice of axes responds to the fact that in the mean-field model σ Y depends inversely on the mechanical coupling strength α; so, the effect of the initial stress distribution are only comparable for different values of α when they are measured relative to the distance between the instability threshold σ c and the dynamical yield stress σ Y (α). In a way, the ratio s d /(σ c − σ Y ) compares the relaxation level (or stress spread) prior to a creep test to the spread for a system flowing very slowly, since σ c − σ Y characterizes the spread of the stress distribution in the zero strain rate limit. When s d and σ c − σ Y become comparable, we would expect σ s Y − σ Y to approach zero. In other words, if the initial aging is not enough (large values of s d ), the overshoot in the stress-strain curve (that distinguishes σ s Y and σ Y ) should cease to be observed. This is consistent with recent observations in particle based simula- tions [51] and confirmed by our data. Both models show . The insets of Figs. 8 and 9 present the creep exponent β s against the relative relaxation. We see, for both models, an increase in β s with increasing level of relative relaxation. Apart from some numerical fluctuations, the collapse of β s obtained for different values of α, shown in the inset of Fig. 9, suggests a master relation between β s and s d /(σ c − σ Y ). The value of β s found at large s d is comparable with experimental measurements [6]. We also plot β against the relative relaxation level s d /(σ c − σ Y ) in Fig.5(b). Curves for different values of α collapse together, while the curve of the elasto-plastic model is a bit shifted aside.

The static yield stress
Previously we have mentioned the relation between the static yield stress σ s Y and the zero strain rate limit of the stress overshoot, we now address this point. We switch to the strain rate controlled protocol and perform shear start-up simulations at different values of strain rate from the same initial conditions used for the creep tests. We record the largest stresses reached in the stress-strain curves produced under different values of strain rates and initial relaxation levels (s d ). We find that for a given s d the stress overshoot decreases with applied strain rate, converging to a finite value when the strain rate goes to zero (see the supplemental material of Ref. [20]). The corresponding limit values for each s d are plotted with crosses in Figs. 8 and 9. We observe that, although the static yield stress and the stress overshoot do not perfectly match, they stay representative each other along the trend as function of s d . Comparing the spatial and the mean-field models, these quantities have a better agreement for small values of s d , which justifies the idea that the underlying physics of the static yield stress observed in creep experiments and that of the quasi-static shear stress overshoot are closely related. We also notice that at large values of s d /(σ c − σ Y ), a small but systematic deviation exists. In such situations, the initial conditions correspond to poorly relaxed systems, with a very short fluidization time, for which deviations from the scenario outlined above may be expected.

Correlations and cooperativity in the creep dynamics
The effective dynamics described by Eqs.9 and 10 and the results shown in Fig. 6 explain, from a mean-field perspective, in terms of populations above and below characteristic values, the underlying physical process involved in creep and fluidization. In the spatially resolved model we are able to further see the spatial distribution of these populations. Therefore, we discuss now the underlying physics of creep dynamics from the point of view of the cooperativity of local plastic events in the elastoplastic model. Let us recall that on each site, the state variable n ij (t) alternates between zero and one respectively for elastic and plastic states. In the creep simulation, one sub-volume or block contributes to the macroscopic strain rate only when in its plastic state so that the time evolution of the state variable field can be used to infer the physical process underlying creep.
In order to implement this, we choose to accumulate the state variable during a given time window The field f ij integrates the plastic activation information during that time window. Here ∆t is chosen to be of the same order of magnitude as t for the plastic activations at a given time scale t to be well represented by f ij (t).
(Actually we have ∆t ≈ 2t/3.) We then compute the following correlation map where · represents ensemble averages. To characterize a "cooperativity level" of the plastic events taking place within a corresponding time window, we define the correlation intensity I C as the integral of the absolute value of the correlation map: This quantity is indicative of how strongly plastic events are correlated during a specific stage of the creep test. The correlation intensity as function of time and the correspondingγ(t) for different values of the imposed step stress are presented in Fig. 10. After a linear increase of both I c (t) andγ(t) for t t mic 1 (not-shown), leading all curves to a common I o c level, we observe distinctive features according to the value of the stress. For the smallest applied stress, the correlation intensity decreases with time and drops to zero, similar to the strain rate. For stresses that lead to steady state flow in the long term, the correlation intensity continues to increase beyond I o c . Note that while the behaviour of the strain rate at early stages is not a clear indicator of fluidisation, an increasing I c (t) appears to be clearly correlated with it.
For those curves leading to a steady flowing state in Fig. 10(a), the corresponding correlation intensities I c (t) (Fig. 10(b)) exhibit two characteristic time scales: (1) a bump of I c before entering the steady state flow, for example between points 3 and 6 on the blue curve and between points 2 and 4 on the red curve. This bump lies roughly in the same range as the inflection point inγ, i.e. corresponds to the fluidization time τ f . This relates the fluidization time scale τ f to a maximum cooperativity of plastic events. The amplitude of the maximum in I c decreases with increasing the applied stress beyond σ s Y . When the applied stress is big enough, less cooperation is needed in the system in order to overcome the static yield stress. The large amplitude in I c observed between points 3 and 6 for the blue curve corresponds to a stress that is very close to the static threshold σ s Y . There, significant spatial correlations are needed in order to fluidize the system. (2) At longer times, after entering the steady state where no other time scales can be recognized from theγ(t) curves, further information comes from the correlation intensity. I c (t) displays a local maximum at long time (point 7 on the red curve) after which it decreases again to a low value in the steady state, point 9. Although the strain rate seems to reach a steady state earlier, the true steady flowing state is achieved only after the correlation level (cooperativity) is relaxed to a steady value. In fact, the strain rate is still evolving up to that moment, but on a very small scale. One recognizes also in Fig. 10 that the steady state cooperativity level of plastic events is inversely related to the steady state strain rate. Intuitively, a fast strain rate tends to activate plastic events rather randomly, while a slow strain rate leaves enough room (time) for correlations in plastic activity to develop [52].

Activity maps and spatial correlation maps
The previous discussion of the development and correlations of cooperativity can be explicitly visualized in our spatial model by snapshots of accumulated plasticity and spatial correlations within a time window, at different stages of the dynamics and for different applied stresses. Figs. 11 shows these maps at the instants marked by the number labels in Fig. 10.
At the beginning, not far from I c (t) = I o c , even when there is significantly more plastic activity for the larger applied stresses than for the smaller one, their correlation maps are almost structureless. A clearer spatial structure of correlation only appears at the onset of fluidization, close to the peak of I c (t) (Fig. 11(b,d)). At this point, the accumulated activity looks homogeneous in space for the large applied stress (Fig. 11(c)), while the cooperativity of plastic events is confirmed by the quadrupolar form in the correlation map displayed in Fig. 11(d). On the other hand, for the small applied stress we already observe the plastic activity organized in several vertical and horizontal bands ( Fig. 11(a)) and a very pronounced vertical correlation pattern in Fig. 11(b), corresponding to the fluidization phase (4 to 5 in Fig. 10). This suggests that a strong correlation developing along the y-direction could be responsible for the burst of plastic events leading to fluidization (and the speed-up ofγ). Of course, there is no a priori preferred direction for cooperativity and the vertical correlation band observed in this example should be equally frequent as an horizontal one. Note that the elasto-plastic model may over-emphasize the stripes formation for the cases of small applied stresses due to the Fourier space implementation of the interaction kernel G and the fact that the model assumes homogeneous elasticity mediating the interactions among plastic events. Nonetheless, transient shear banding during creep is indeed observed in more realistic MD models [4].
After the fluidization regime, the behavior of plastic activity for the smaller applied stress is similar to that for the larger applied stress (not shown). For the large applied stress, the quadrupolar correlation gradually develops until one arrives at the last maximum in I c and drops back to a very weak spatial pattern in the steady state. Correspondingly, plastic activity is more organized into thin slip lines at the last peak of the correlation intensity than in the steady state. The phenomenology is a bit more complex for the small applied stress, close to σ s Y . The burst of cooperative plastic events that organize into the shear bands in the fluidization regime not only speeds up the strain rate but may also induce sites outside the bands to become unstable: since plastic activations inside a band tend to decrease the stress, this should be balanced by an increase of local stresses outside the bands in order for the overall averaged stress to keep constant. As a result, more plastic events take place randomly outside the shear bands; this acts against overall cooperativity: the correlation intensity drops significantly after the large bump (3 to 6 in Fig. 10). The remaining process afterwards is very much like the one in the case of the large applied stress. Activity maps show thin random slip lines for both stresses and the correlations are both very weakly quadrupolar. The only difference is that the last maximum in I c (t) appears later and the final stationary correlation pattern is a bit more pronounced for the smaller stress.

CONCLUSIONS
We used both spatially-resolved and mean-field mesoscopic models to study the creep behavior of athermal amorphous materials with different initial relaxation degrees. Despite the simplicity of the models, they are sufficient to reproduce the S-shaped strain rate response observed in experiments. Further, the two models are consistent in both qualitative and quantitative manners, so that the mean-field model can be considered as a simplified way to understand the more realistic spatially resolved model. We measured the power law slowing dowṅ γ ∼ t −µ in the creep regime. We found that the exponent µ produced by both models lies in the same numerical range as experimental results but is not universal with respect to the applied stress and the level of initial relaxation. We distinguished, within the framework of the models, the different underlying physical processes for the two time scales τ m and τ f that characterize the creep behavior. τ m is determined by the initial stress distribution P 0 (σ) around the marginal stability threshold σ c , while τ f is closely related to subsequent plastic activations and spatial cooperativity of plastic events. We interpreted our results on the relation between the fluidization time τ f and the applied stress σ EXT by defining a static yield stress σ s Y , which increases with initial relaxation. A convincing power law τ f ∼ (σ EXT − σ s Y ) −βs is observed, with β s increasing when shortening the initial aging and taking values comparable to those reported in experimental studies. Finally, we defined an intensity of spatial cooperativity that can serve as a precursor to distinguish systems that fluidize from those stuck at the early stage of the creep phase. The onset of the fluidization regime is associated, especially for small stresses, with a strong spatial cooperativity. Moreover, we noticed that spatial correlations in cooperativity seem to be qualitatively different between systems that undergo a creep regime prior to the fluidization and those that fluidize directly.