Statistical dynamics of early creep stages in disordered materials

TL;DR

This paper investigates the early stages of creep deformation in disordered materials using a mesoscopic elastoplastic model, revealing how structural evolution influences event statistics and exhibits Omori law behavior.

Contribution

It introduces a stochastic elastoplastic model to analyze structural and statistical evolution during creep, linking microstructural heterogeneity to observed laws.

Findings

Event rate follows Omori law, related to Andrade's creep law.

Transition from Poisson to power-law inter-event times during creep.

Microstructural heterogeneity determines the exponents of observed laws.

Abstract

When materials are loaded below their short-term strength over extended periods, a slow time-dependent process known as creep deformation takes place. During creep deformation, the structural properties of a material evolve as a function of time. By means of a generic coarse-grained mesoscopic elastoplastic model which envisages deformation as a sequence of stochastically activated discrete events, we study the creep deformation of disordered materials. We find that the structural evolution of the material during creep modifies not only the average material properties but also changes the statistics of those properties. We analyze the emergence of correlations in the strain localization and deformation activity patterns, the variation of the event rate and the evolution of the inter-event time distribution. We find that the event rate follows the Omori law of aftershocks, which is the discrete counterpart of Andrade's transient creep law, and that the exponent of these laws only depends on the microstructural heterogeneity. Finally, we find during the initial stages of transient creep a transition from Poisson distributed inter-event times towards a non-trivial power law distribution.

Equations

$$\Phi =\sqrt{(3/2)\mathbf{\Sigma }{}^{\prime }:\mathbf{\Sigma }{}^{\prime }}-\hat{\Sigma }={\Sigma }_{\mathrm{eq}}-\hat{\Sigma }.$$

$$\Delta \mathbf{ϵ}=\widehat{\mathbf{ϵ}}\Delta ϵ$$

$${\widehat{ϵ}}_{ij}=\partial \Phi /\partial {\Sigma }_{ij}.$$

$$\nu (r_i)={\nu }_{\mathrm{el}}\exp \left(-\frac{E(\mathbf{\Sigma })}{k_{\mathrm{B}}T}\right).$$

$$\nu (r_i)={\nu }_{\mathrm{el}}\exp \left(\frac{\Phi (r_i)}{{\Sigma }_T}\right)$$

$$P(\hat{\Sigma }|\lambda ,k)=\frac{k}{\lambda }{\left(\frac{\hat{\Sigma }}{\lambda }\right)}^{k-1}\text{exp}\left(-{\left(\frac{\hat{\Sigma }}{\lambda }\right)}^k\right)$$

$$\rho (ϵ)=\frac{⟨x{x}^{\prime }{⟩}_{ϵ}-⟨x{⟩}_{ϵ}^2}{\sigma (x,ϵ{)}^2}$$

Full text

††preprint: APS/123-QED

I Introduction

It is well-known M.Sc. (1905); Andrade and Trouton (1910) that when materials are loaded below their short-term strength over extended periods, a slow time-dependent deformation known as creep takes place. Understanding creep remains of general concern in view of the long-term behavior of structural materials Lee et al. (1991) and geosystems such as Earth crust faults zones Sleep and Blanpied (1992). Creep deformation is characterized by several deformation regimes. Upon application of an external load, an instantaneous deformation takes place which is followed by a regime of decelerating deformation rates known as transient or primary creep. Transient creep is described by the Andrade law Andrade and Trouton (1910); Louchet and Duval (2009), which relates deformation to the time since application of the external load as $\epsilon\sim t^{m}$, with $m\approx 2/3$. Nonetheless, the value of $m$ is not universal and depends on the specific material under study Louchet and Duval (2009). Alternatively, transient creep is sometimes reported to be logarithmic, $\epsilon\sim\textrm{log}(t)$ M.Sc. (1905); Cottrell (1997). After the transient regime, if the temperature is high enough, deformation enters the stationary or linear creep regime, $\epsilon\sim t$ Louchet and Duval (2009); Main (2000).

The description of transient creep in terms of the Andrade law has traditionally been associated with a conceptualization of creep in the framework of continuum mechanics, where plastic deformation appears as a smooth and deterministic process. Nonetheless, over the last decades, it has become clear that plasticity on small scales in which individual plastic events can be resolved appears as a stochastic and intermittent phenomenon (Zaiser, 2013; Schneider et al., 2009). When individual deformation events can be resolved, the slowing-down of the intermittent activity is described by the Omori law. The Omori law was originally established for earthquakes Omori (1894); Utsu et al. (1995) even before the early observations of transient creep which led to the Andrade law Andrade and Trouton (1910). Specifically, it describes the rate of earthquake aftershocks as $\dot{n}\sim t^{-p}$, where $\dot{n}$ represents the rate of events and $p$ is typically close to 1. Omori-like behavior of plastic deformamtion activity indicators has been observed in a wide diversity of materials such as, e.g., rocks Lennartz-Sassinek et al. (2014), porous materials Baró et al. (2013), protein gels Leocmach et al. (2014) or bulk metallic glasses McFaul et al. (2018).

Given the underlying fluctuating nature of the plastic activity, relating the phenomenological laws of creep to specific microstructural processes has been a major task in understanding creep deformation. In the case of crystals, this has been achieved by means of dislocation theory Cottrell (2004); Miguel et al. (2002), where the elementary deformation events correspond to the activation of dislocation segments and their motion between metastable configurations, and the structural disorder arises from the complex stochasit patterns that are formed by the dislocation system Zaiser and Hähne (1999); Zaiser and Aifantis (2006). In the case of metallic glasses, the elementary deformation events correspond to shear transformations Argon (1979); Schuh et al. (2007). A shear transformation can be understood as a local atomic rearrangement encompassing around 100 atoms which re-arrange in a shear-like fashion in order to accommodate the local shear stress. The structural disorder inherent to a metallic glass is the result of non-equilibrium features frozen in the structure after a liquid melt is quenched into the glassy state.

Similar creep behavior is found across a wide range of materials, such as polycrystalline metals Andrade and Trouton (1910), metallic glasses Krisponeit et al. (2014); McFaul et al. (2018), rocks Heap et al. (2011), paper Rosti et al. (2010) or granular matter Nguyen et al. (2011). Such independence of microstructural details suggests a conceptualization of creep phenomena within a framework that does not rely on material-specific microstructural processes. To this end, mesoscale elastoplastic models of disordered materials have been widely employed in recent years to reproduce and analyze intermittent plastic activity from a coarse-grained perspective Nicolas et al. (2018). In such models, the material is discretized into a grid of mesoscopic elements of size equal to or above the length scale of the elementary deformation events, with statistically distributed elemental properties in order to represent the underlying disordered microstructure.

In our work, we consider a mesoscale elastoplastic model based on the concept of local yield threshold, which corresponds to a mesoscopic internal state variable akin to the yield stress of continuum mechanics and characterizes the local resistance to plastic deformation. Once the locally acting shear stress overcomes the yield threshold of a mesoscopic element, a localized plastic event is triggered within that element. This process aims at capturing the essentials of plastic flow at the microscale without resolving atomistic details. The concept of yield threshold is of broad applicability, irrespective of the nature of the underlying microstructure. Advantages of the concept reside in the possibility of directly measuring the yield stress probability distribution by probing mesoscopic regions in molecular dynamics Barbot et al. (2018); Patinet et al. (2016) or DEM simulations, and by the connection it makes between plasticity and interface depinning problems Budrikis et al. (2017); Lin et al. (2014); Ozawa et al. (2018); Talamali et al. (2011).

Although mesoscale elastoplastic models have successfully reproduced much of the phenomenology associated with plastic deformation of disordered materials under strain or stress driven conditions Budrikis et al. (2017); Budrikis and Zapperi (2013); Sandfeld et al. (2015); Talamali et al. (2011); Lin et al. (2014); Tüszes et al. (2017), comparatively little attention has been given to creep deformation Merabia and Detcheverry (2016); Castellanos and Zaiser (2018); Bouttes and Vandembroucq (2013); Liu et al. (2018). In this paper, we implement a model that is closely related to the model of creep failure by Castellanos and Zaiser (2018), but focus on the short-time transient dynamics and the long-time stationary flow regime. The model introduces stochastic behavior by considering thermally activated events and statistically distributed microstructural properties, which allows us to link the changes in the macroscopic creep response to the structural evolution. The model is able to capture the aging of the structure as well as mechanical rejuvenation as a function of time, external stress, temperature, and microstructural disorder. We study the spatial correlations between plastic events, strain localization, event rate and inter-event waiting time distributions during transient and stationary creep regimes.

II Model

We coarse grain the microscopic details of plastic deformation events and represent the material as a 2D array of yielding elements Budrikis et al. (2017); Budrikis and Zapperi (2013); Sandfeld et al. (2015); Talamali et al. (2011); Lin et al. (2014); Tüszes et al. (2017) of linear size $L$ and volume $V_{\rm el}$ centered at positions ${\bm{\mathnormal{r}}}_{i}$. The state of an element $i$ is represented by (i) a local stress tensor $\bm{\Sigma}({\bm{\mathnormal{r}}}_{i})$ which is the superposition of the stresses resulting from external boundary conditions (for creep: temporally constant applied tractions) and internal stresses resulting from the heterogeneity of the plastic strain field, (ii) its accumulated plastic strain $\bm{\epsilon}({\bm{\mathnormal{r}}}_{i})$ and (iii) a local yield threshold $\hat{\Sigma}({\bm{\mathnormal{r}}}_{i})$ which characterizes the internal state of an element. Plastic deformation is assumed to be governed by the yield function

[TABLE]

Here, $\Sigma{{}^{\prime}}$ is the deviatoric part of the stress tensor and the equivalent stress $\Sigma_{\rm eq}$ is defined in such a manner that the model amounts to a stochastic generalization of J2 plasticity (for further generalizations see Budrikis et al. (2017)). The yield function depends on the local yield threshold and the local stress tensor, which is computed from the external boundary conditions and the plastic strain field $\bm{\epsilon}({\bm{\mathnormal{r}}}_{i})$ using standard Finite Element methodology.

On microscopic scales below the scale of resolution of our model, microscopic plastic re-arrangements are assumed to result in deformation events which produce on the element scale a tensorial plastic strain increment

[TABLE]

where $\Delta\epsilon$ is chosen such that the local equivalent stress is reduced by a factor $0<e\leq 1$. The tensor $\bm{\hat{\epsilon}}$, which gives the ’direction’ of the local strain increment is in the spirit of an associated flow rule chosen to maximize energy dissipation by setting

[TABLE]

II.1 Rules for local activation of deformation events

The plastic flow law we adopt is characterized by the idea that deformation occurs as a stochastic sequence of local events. These events are activated according to the local values $\Phi({\bm{\mathnormal{r}}}_{i})$ of the yield function. We use the following rules:

• •

An event is activated instantaneously if $\Phi({\bm{\mathnormal{r}}}_{i})>0$. We denote this process as mechanical activation.

• •

If $\Phi({\bm{\mathnormal{r}}}_{i})<0$, an event is activated with finite rate that depends on temperature $T$ according to

[TABLE]

where $\nu_{\rm el}$ is an attempt frequency for event activation within the element volume $V_{\rm el}$. We approximate the stress dependence of the characteristic activation energy $E$ by a linear dependency on the equivalent stress, $E=E_{0}-V_{\rm A}\Sigma_{\rm eq}$where $V_{\rm A}$ is an activation volume. The activation barrier goes to zero if $\Phi=0$, hence $E_{0}$ relates to the activation threshold stress via $E_{0}=\hat{\Sigma}V_{\rm A}$ and we can write the activation rate alternatively as

[TABLE]

where the parameter $\Sigma_{T}=k_{\rm B}T/V_{\rm A}$ characterizes the influence of thermal fluctuations on event activation. Besides thermal activation over a barrier, this parameter might also stem from other thermally activated processes that may trigger an event, e.g. the rate of chemical attack in chemically assisted microcracking. Beyond thermal activation, the same parameter might also be interpreted in terms of an effective magnitude of stochastic stress fluctuations of non-thermal origin affecting the region of interest, as long as such fluctuations can be represented by an effective temperature. Irrespective of the physical origin of the fluctuations, we denote this activation process as thermal activation.

• •

The duration of a deformation event is assumed negligibly small.

In the limit where the activation thresholds $\hat{\Sigma}$ are spatially uniform and the stress changes $\Delta\epsilon$ are infinitesimally small, our model reduces for $T\to 0$ (no thermal effects) to a standard ${\rm J}_{2}$ plasticity model. On the other hand, at low external stress, the model reduces to a viscoplastic creep model where the rate of plastic flow is given by $\dot{\bm{\epsilon}}=\bm{\hat{\epsilon}}\nu_{0}\Delta\epsilon\exp(\Phi/\Sigma_{T})$.

II.2 Distribution of yield thresholds

Statistical heterogeneity of the material is represented by considering the local yield thresholds $\hat{\Sigma}({\bm{\mathnormal{r}}}_{i})$ as random variables. Assuming the weakest-link hypothesis at the microscopic scale Alava et al. (2009), a mesoscopic element yields when the weakest of the microscopic elements it coarse-grains yields. In this case the mesoscopic yield thresholds are expected to follow a Weibull distribution. We consider thus a Weibull probability density with a mean value of $\hat{\Sigma}_{0}$ and exponent $k$

[TABLE]

where the scale parameter $\lambda$ is chosen to ensure that the distribution has a mean of $\hat{\Sigma}_{0}$ for a specific given value $k$. The value of $k$ relates to the width of the threshold distribution and consequently characterizes the disorder of the yield thresholds. Specifically, lower values of $k$ correspond to broader distributions, i.e. higher disorder. After each deformation event, the local yield threshold of the deforming element is assigned a new value from the probability distribution (6) in order to represent structural changes. For a generalization involving strain softening, see Castellanos and Zaiser (2018).

II.3 Simulation protocol

Simulations are carried out under pure shear conditions by imposing on the free surfaces of the system spatially uniform tractions giving rise to homogeneous shear stress, in the following denoted as $\Sigma$, which is kept fixed during a simulation. Deformation proceeds as a sequence of plastic deformation avalanches. A thermally activated event initiates an avalanche in form of a cascade of mechanically triggered events happening on a short time scale which we assume to be negligible in comparison with the inter-event time between different thermally activated events. We avoid thus competition between time scales Liu et al. (2016), which is beyond the scope of the present investigation. The thermal events are selected by the Kinetic Monte Carlo Method with the local activation rates given by Eq. (5); the KMC formalism yields both the location of a the thermally activated event and the inter-event time that has elapsed since the last avalanche. After event initiation, we increase the local strain at the activated site according to Eqs. (2,3). We then re-compute the stress field and check whether, as a consequence of stress re-distribution, the condition $\Phi({\bm{\mathnormal{r}}}_{i})<0$ is fulfilled at any site. These sites also become activated and yield, leading to further stress changes and possible activation of further elements. The ensuing avalanche proceeds adiabatically in a series of deformation steps in each of which one or more elements are activated and yield (parallel update) until the inequality $\Phi({\bm{\mathnormal{r}}}_{i})<0$ is not fulfilled anywhere Budrikis and Zapperi (2013); Budrikis et al. (2017). The avalanche size $S$ is then defined as the total number of mechanically activated events between avalanche initiation and termination. After termination of the avalanche, a next Kinetic Monte Carlo step determines the initiation site and initiation time of the next avalanche. We repeat this process until a specified end strain is reached. To obtain statistically representative results, we perform ensemble averages over many realizations of the disorder (typically $10^{3}$ samples for each set of parameters).

II.4 Parameters

In the simulations, stress is measured in units of $\hat{\Sigma}_{0}$, the strain in units of $\hat{\Sigma}_{0}/E$ ( where $E$ is Young’s modulus) and time in units of $\nu_{0}^{-1}$. The model depends on the non-dimensional coupling constant $1-e$ which controls the characteristic fraction of the local stress that is elastically re-distributed after an event. The model relates to particular disordered systems through the values of $e$, $\hat{\Sigma}_{0}$ and $k$. Unless otherwise stated, during this work we assume the default values $e=0.95$, $L=64$, $\Sigma_{T}=0.0075$ and $k=4$.

Since creep loading requires the application of a constant external load, we need to determine the value of such load before performing creep tests. To this end, we use as reference the critical macroscopic yield stress $\Sigma_{\rm c}$ at which the system enters a state of sustained macroscopic plastic flow without the need for thermal activation. The value of $\Sigma_{\rm c}$ averaged over many simulations for different system parameters is shown in Fig. 1. Specifically, for the default set of parameters given above, we find $\Sigma_{\rm c}=0.668$. We set a default value for the external stress of $\Sigma=0.7\Sigma_{\rm c}$ during our creep simulations unless otherwise stated.

III Simulation results

III.1 Creep curves

Based on the evolution of the average strain as a function of time (see Fig. 2, top), we observe two deformation regimes. The first regime, referred to as transient creep, is characterized by a decelerating strain rate. The second regime, referred to as stationary creep, is characterized by a constant strain rate. Tertiary creep, which is characterized by an accelerating strain rate prior to failure, is not in the focus of the present work but has been discussed in detail in a parallel study (Castellanos and Zaiser, 2018).

Immediately after the application of the load at $t=0$, all elements with a yield threshold lower than the external stress yield in an instantaneous and uncorrelated manner. Consequently, a macroscopic avalanche with a characteristic size of order $\langle S\rangle\sim L^{2}$ takes place. This avalanche can be interpreted as the instantaneous creep stage Louchet and Duval (2009). After the instantaneous creep, the average avalanche size quickly drops as activity is mainly due to thermal activation of individual plastic events. This gives rise to the transient regime, during which deformation continues to decelerate. After the transient regime, the deformation enters a stationary regime, where activity is controlled by thermally activated events, with an amount of mechanically triggered sequels that depends on the stress level and diverges as $\Sigma\to\Sigma_{c}$ but is typically small away from this critical point.

In the following, we consider different magnitudes that will help us gain insight into the internal dynamics of systems undergoing creep deformation. We study the evolution of such magnitudes along the creep process and their dependence on system parameters and deformation conditions.

III.2 Evolution of yield thresholds

We first consider the time evolution of the average local yield threshold $\langle\hat{\Sigma}\rangle$, and the threshold scatter measured by the standard deviation of thresholds ${\rm std}(\hat{\Sigma})$. As shown in Fig. 3 (top) and Fig. 4 (top), the average threshold first increase from its initial value of $\langle\hat{\Sigma}\rangle=1$. This hardening is responsible for the slowing down of creep deformation in the transient regime. Commonly known as statistical hardening, it is the consequence of a survival bias due to the statistically distributed yield thresholds: The lowest threshold elements are, according to the rules of Sec. II.1, more likely to yield and the re-assignment of a new threshold leads, on average, to a higher local threshold. The evolution of the average threshold can be interpreted as structural aging, by which thermally activated plastic activity leads to a more stable structure. The value of ${\rm std}(\hat{\Sigma})$ is found to drop until reaching a minimum value, after which it grows until attaining its plateau value at the stationary creep regime (see Fig. 3 and see Fig. 4, bottom). The existence and location of the minimum depends on temperature (it disappears at high temperatures) but is robust upon variation of the external stress.

Statistical hardening saturates as the thresholds move to the high strength side of the distribution (6) until a dynamic equilibrium situation is reached where the mean and scatter of the yield thresholds in the sample no longer change. Thus, the system enters the stationary creep regime which is characterized by a constant deformation rate, as shown in Fig. 2 (bottom). We denote the value of the stationary average yield threshold as $\langle\hat{\Sigma}_{\rm II}\rangle$ and that of the threshold scatter as ${\rm std}(\hat{\Sigma}_{\rm II})$. Fig. 5 shows the values of $\langle\hat{\Sigma}_{\rm II}\rangle$ for different system parameters, and Fig. 6 shows the coefficient of variation ${\rm std}(\hat{\Sigma}_{\rm II})/\langle\hat{\Sigma}_{\rm II}\rangle$. We find that, as the external stress is increased, the average threshold stabilizes at a lower value while the coefficient of variation increases. Therefore, an increase of the external stress reverts to some extent, the effects of aging. This stress dependence of structural stability is commonly known as mechanical rejuvenation Gruber (1978); Zhang et al. (2017); Warren and Rottler (2008). Increasing temperature is found to enhance the effect of mechanical rejuvenation. We find that systems with higher initial disorder (i.e., lower $k$) reach higher saturation values of the local strength. This is in line with the observations of Tüszes et al. (2017) who report, for systems loaded in displacement control, that strength increases with increasing degree of disorder. The system size is found to have no discernable influence on the statistical behavior.

III.3 Spatial activity patterns

III.3.1 Strain localization

During the initial stages of the transient regime, internal stresses are almost absent, and the system behavior is controlled by externally applied shear stress which is spatially homogeneous. As a result, plastic deformation is activated in a spatially uniform random manner (see Fig. 2, patterns). As deformation proceeds, the internal stress field created by previous plastic events induces elastic couplings within the system, which leads to self-organization of the deformation activity and the emergence of strain patterns in the form of non-permanent, mutually perpendicular shear bands. The orientation of the bands is the result of the loading mode: our external load induces a pure shear state with principal axes oriented along $\pm\pi/4$ with respect to the horizontal axis in Figs. (2) and (8). The plastic events follow this orientation, and, as a consequence, positive stress redistribution favoring the formation of bands takes places along the directions [math] and $\pi/2$ with respect to the horizontal axis. However, the ensuing strain localization is only transient: As the creep strain increases during the stationary creep regime, the plastic strain pattern arising from the superposition of strain bands remains statistially homogeneous.

III.3.2 Event correlation

To quantitatively study the spatial activity, we consider the Pearson correlation coefficient $\rho$ between the locations of the triggering events of two consecutive avalanches. Specifically, we consider the projected locations over the horizontal axis, $x$, and $x^{\prime}$. As discussed by Castellanos and Zaiser (2018), care must be taken in the case of the symmetry breaking induced by, e.g., permanent strain localization in the form of a macroscopic shear band. However, since we study the stationary creep regime without the presence of any softening mechanism, the activity is symmetric with respect to the projection over any pair of perpendicular directions. We defined thus the coefficient

[TABLE]

where $\langle\cdot\rangle_{\epsilon}$ denotes an average over a narrow strain window centered at $\epsilon$ and $\sigma(\cdot,\epsilon)$ is the standard deviation of event locations within that strain window (when no $\epsilon$ is present, the operation is performed over all the complete simulation). The evolution of the correlation coefficient $\rho$ with strain is shown in Fig. 2 (bottom). The initial lack of elastic coupling is reflected by near-zero correlation coefficient $\rho$, while on the approach to the stationary regime the coefficient saturates at a value $\rho_{\rm II}$.

The dependence of $\rho_{\rm II}$ on different system parameters is shown in Fig. 7. Since spatial correlations are the result of the internal stress field induced by the plastic events, a higher (spatially homogeneous) external stress reduces the effects of such internal stress field and leads to a drop of spatial correlation. Similarly, by increasing temperature we promote the stochastic yielding of elements, leading to reduced spatial correlation. The drop as the system size increases stems from the fact that consecutive events are more likely to occur far away from each other in bigger systems. We observe that, as the disorder increases, the correlation increases. This behavior can be understood by looking at the strain patterns obtained at the same value of strain with different values of $k$. As shown in Fig. 8, a higher disorder (lower $k$) promotes activity localization in small clusters that are uniformly scattered across the system. In the case of a lower disorder, plastic activity gives rise to comparatively longer-wavelength features, which leads to a lower correlation coefficient.

III.4 Event rate

The event rate $\dot{n}$ is defined as the number of avalanches per unit time. Since avalanches are considered as instantaneous events and are triggered by thermally activated events, $\dot{n}$ matches the rate of thermal activation.

III.4.1 Transient regime

The event rate is found to decrease after the application of the external load as $\dot{n}\sim t^{-p}$ with $p\approx 0.89$ (see Fig. 9, left). This relation is known as Omori law and describes the decay of activity with time after a large perturbation occurs in the system. It is commonly used in geophysics to model the sequence of aftershocks after a main shock (Utsu et al., 1995; Schmid and Grasso, 2012) and has been found in materials failure tests after a partial rupture of the sample (Baró et al., 2013). In our case, the perturbation corresponds to the sudden application of the external load at $t=0$. The drop of event rate can be explained as the result of the growth of the yield thresholds (Fig. 2, bottom) which leads to lower activation rates in Eq. 5.

We have studied the parameter dependence of the Omori exponent $p$. We find that, within error bars, the value of $p$ does not vary with the value of the applied external stress, the temperature or the system size. However, we find a systematic evolution of the value of $p$ upon variations of the Weibull exponent $k$ (see Fig. 9, right). The exponent $p$ decreases upon increasing the threshold disorder (i.e., lower $k$), and seems to converge towards $p=1$ for $k=0$. Since avalanches have, during our simulations, a constant characteristic size $\langle S\rangle$, the cumulative deformation during the creep process can be obtained by integrating with respect to time the event rate, $\epsilon\sim\int\dot{n}dt$. Consequently, we recover the Andrade law $\epsilon\sim t^{m}$, with an exponent $m=1-p$, in the case of highly ordered systems of low $k$. On the other hand, as disorder increases and $p\to 1$, we recover logarithmic creep $\epsilon\sim\textrm{log}(t)$.

III.4.2 Stationary regime

When the creep process enters the stationary regime, the average threshold becomes constant (Sec. III.2). Since the element activation rates are related to the thresholds through Eq. 5 and Eq. 1, the activity in the stationary regime $\dot{n}_{\rm II}$ becomes eventually constant. The impact of different simulation parameters on $\dot{n}_{\rm II}$ is shown in Fig. 10. The exponential dependence on applied stress and inverse temperature are direct results of the underlying thermal activation process as given by Eq. 5, i.e., our model is in accordance with viscoplastic creep models. The event rate scales with system size as $L^{2}$, which is the consequence of considering thermal activation a Poisson process in which $L^{2}$ elements with the same average activation rate (in the stationary regime) attempt to yield simultaneously. Increasing the disorder is found to strongly reduce the event rate in the stationary regime, which is explained by the the higher saturation values of the local strength reached (Fig. 5) and again illustrates the paradigm ’more disordered is stronger’.

III.5 Inter-event time distributions

The inter-event time between two thermally activated events is computed by the Kinetic Monte Carlo method. Thus, thermal activation is envisaged as a Poisson process in which the mesoscopic elements attempt to yield independently. However, the distribution of activation barriers is in itself an emergent property as the local activation barriers $\Phi$ depend on the internal stresses and therefore change due to the mutual interactions between events occurring at difffernt times. Therefore, the probability density $P(\Delta t)$ is expected to reflect the different degrees of activity correlation found along the creep curve.

Fig. 11 shows the evolution of $P(\Delta t)$ during the initial stages of the transient regime. Initially, the probability distribution is of exponential shape, suggesting uncorrelated activity. This is in line with the initial lack of spatial correlation of plastic activity of Sec. III.3.2. As deformation proceeds, the average threshold grows (Fig. 3), and the system becomes thus more stable. Consequently, the average time for a thermal fluctuation to occur which can overcome an energy barrier increases, which is reflected in the shift of the distribution cut-off toward bigger values. Still during the transient regime, the probability distribution becomes dominated by a power-law regime with an exponent close to $-1$, incompatible with Poissonian statistics. This happens in parallel with the growth of spatial correlations after the application of the load (Fig. 2, bottom). During the stationary creep regime, the distribution of inter-event times becomes stationary and retains its power law charateristics.

We note that, as shown by Castellanos and Zaiser (2018), during the third creep stage, i.e. during the approach to failure, the inter-event time distribution undergoes exactly the reverse evolution - it changes from a power law with exponent close to -1 into an exponential distribution, and this evolution is concomitant with a loss of spatial correlation between subsequent events.

IV Discussion and Conclusions

As described in Sec. III.2, the transient creep regime is associated with hardening, reflected in the growth of the average yield threshold, which in turn is the result of the exhaustion of low-threshold elements. The saturation value of the average yield threshold in the stationary creep regime is found to depend on the system parameters (Fig. 5). At low temperatures, low threshold elements are typically chosen for yielding. The renewed threshold is thus likely to be significantly higher than the previous one, which enhances statistical hardening. On the other hand, a higher temperature allows for more stable elements to be chosen for yielding, reducing thus the effects of statistical hardening. The stationary threshold is determined, at a certain temperature, by the dynamical equilibrium achived between both processes, namely the seletion of low-threshold elements which become stronger after yielding and the selection of strong elements which become softer after yielding. Increasing the external stress has an effect similar to that of increasing temperature. This can be understood by noticing that the stress lowers the activation barriers in a global fashion. Consequently, at a higher external stress, strong elements are chosen for yielding more often than they would at low external stress. This has the effect of shifting the dynamical equilibrium of thresholds toward lower values. By increasing the threshold disorder (i.e., reducing $k$), the stationary average yield threshold is increased, which is a direct consequence of fattening the high-strength tail of the threshold probability distribution.

We find that the decrease of activity during the transient creep regime follows an Omori-like law with an exponent $p$ which is generally close to $1$, as is experimentally well-known for a wide range of materials, and recently established in the specific case of metallic glasses McFaul et al. (2018). Specifically, we observe that $p$ does not depend on the external stress, temperature or system size but varies systematically with the degree of disorder. The exponent seems to converge towards $p=1$ at high disorder, while the value decreases with the disorder. As pointed out in Sec. III.4.1, the variation of $p$ with microstructural disorder might explain the variability observed in the exponent of the Andrade law of transient creep or the existence of logarithmic creep. Moreover, the exclusive dependence of the exponent on the disorder of the yield thresholds might allow for the quantitative estimation of the disorder at the mesoscopic scale. The disorder at the mesoscale is a crucial parameter in the modeling of the plastic activity of disordered materials Tüszes et al. (2017); Budrikis et al. (2017); Popović et al. (2018); Ozawa et al. (2018), which is on the other hand difficult to establish experimentally. It is therefore of interest establishing links between the macroscopic observable response and a quantitative measure of structural heterogeneity.

The evolution of the inter-event time distribution as deformation proceeds (Fig. 11) is characterized by a transition from an initially exponential shape to a power-law-like distribution with an exponent of approximately $-1$. This value is in agreement with experimental findings on the creep deformation of bulk metallic glasses Krisponeit et al. (2014); McFaul et al. (2018). The transition from uncorrelated to correlated activity suggested by the distributions of inter-event times is supported by the correlation coefficient of spatial activity (Sec. III.3.2). Interestingly, Krisponeit et al. (2014) found that the distribution of inter-event times evolves with time, from a distribution initially dominated by a cut-off towards a power-law-like distribution with an exponent of approximately $-0.8\pm 0.1$, which agrees well with our observations.

In conclusion, we have presented a mesoscale elastoplastic model which can reproduce the time evolution of the intermittent plastic activity of disordered materials under creep conditions. We find that the structural evolution as a function of time is concomitant to a variation of the spatial correlation of events, to the emergence of strain localization patterns, and to the variation of the event rate and the inter-event times distribution. In summary, our approach allows us to study the creep dynamics from a holistic perspective, establishing contact between diverse approaches which are frequently considered separately, the macroscopic smooth flow response, the microscopic stochastic intermittent response, and the evolution of the structural properties as expressed by the distribution of local activation thresholds.