# Modeling the creep properties of olivine by 2.5-dimensional dislocation dynamics simulations

Francesca Boioli, Philippe Carrez, Patrick Cordier, Benoit Devincre, Matthieu Marquille

► To cite this version:

Francesca Boioli, Philippe Carrez, Patrick Cordier, Benoit Devincre, Matthieu Marquille. Modeling the creep properties of olivine by 2.5-dimensional dislocation dynamics simulations. Physical Review B: Condensed Matter and Materials Physics (1998-2015), 2015, Physical Review B: Condensed Matter and Materials Physics (1998-2015), 92 (1), pp.14115. 10.1103/PhysRevB.92.014115. hal-03390684

HAL Id: hal-03390684

https://hal.univ-lille.fr/hal-03390684

Submitted on 21 Oct 2021

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Modeling the creep properties of olivine by 2.5-dimensional dislocation dynamics simulations

Francesca Boioli,\* Philippe Carrez, and Patrick Cordier

UMET, University of Lille1, Villeneuve d'Ascq, France

Benoit Devincre

LEM, CNRS-ONERA, Chatillon, France

Matthieu Marquille

CRI, University of Lille1, Villeneuve d'Ascq, France

(Received 16 December 2014; revised manuscript received 10 March 2015; published 23 July 2015)

In this work we performed 2.5-dimensional (2.5D) dislocation dynamics simulations coupling climb with the glide dislocation motion to model the creep behavior of olivine, one of the main component of the Earth's upper mantle. In particular, we present an application of this method to determine the creep strain rate in a material with high lattice resistance, such as olivine. We show that by including the climb mechanism we reach steady state creep conditions. Moreover, we find that a creep power law with a stress exponent close to 3 can be extracted from our simulations and we provide a model based on Orowan's law to predict the creep strain rates in the high temperature and low stress regime. The model presented is relevant to describe the plastic flow of olivine in the Earth's mantle deformation conditions and can be useful to derive the high temperature creep behavior of other materials.

DOI: 10.1103/PhysRevB.92.014115

PACS number(s): 61.72.Hh, 61.72.Ff, 62.20.fq, 91.60.Dc

# I. INTRODUCTION

Large scale flow in the Earth's mantle involves plastic deformations of rocks and their constitutive minerals. Due to the extremely slow strain rate conditions in the Earth's mantle (~$10^{-14}$ s⁻¹), it is very challenging to identify the fundamental mechanisms controlling the plastic behavior of its constitutive minerals. Thus, the development of a multiscale approach linking the atomic scale properties and the microscopic elementary mechanisms to the macroscopic behavior is needed [1]. One of the key steps in this approach is the description of dislocation-based intracrystalline plasticity. To this aim, discrete dislocation dynamics (DD) has been developed to describe the collective behavior of dislocations at the mesoscale, filling the gap between atomic scale properties and plastic behavior at the macroscopic scale [2,3]. At low homologous temperatures ($T < 0.4T_m$), DD simulations have proved to give accurate descriptions of many plasticity problems in metals [4,5], semiconductors [6,7], as well as in an oxide, MgO [8], and a silicate, olivine [9]. In this regime, the displacements of dislocations, which are the main carriers of elemental plastic deformation, is mainly controlled by the glide mechanism. On the contrary, at high temperatures ($T > 0.4T_m$), as diffusion processes become important, dislocations can also move via the climb mechanism, which is a nonconservative dislocation motion and requires the absorption or emission of point defects [10,11]. This mechanism is generally not included in classical DD formulations [2,3,12], since most of the previous

DD studies have been focused on the low temperature regime. In the last years however, some DD formulations have included coupled glide and climb mechanisms to model the annealing of dislocation loops [13,14]. Following a similar approach, where the climb velocity is directly related to the flux of point defects absorbed or emitted by dislocations, climb has been incorporated in a two-dimensional (2D) DD code in order to describe the Al mechanical behavior at elevated temperature in micropillars [15], thin films [16] and single crystal [17], and polycrystals [18]. Recently a new 3D formulation of the DD method including climb has been proposed by Po et al. [19]. Although it is generally accepted that climb strongly affects the dislocations microstructure and evolution at high temperature, for example, allowing dislocations to bypass obstacles, or to reduce the back stresses induced by dislocation pileups on dislocation sources, the effect of climb on the creep behavior is still poorly known in many materials.

In this work we focus on the creep behavior of olivine at high temperature and ambient pressure, under low applied stresses. Olivine ($Mg,Fe)_2SiO_4$ is the most abundant and the weakest constituent of the Earth's upper mantle in a wide range of conditions. As a consequence, predicting its mechanical behavior is of fundamental importance to model the rheology of the upper mantle. At high temperatures, both climb and cross slip processes may play an important role on creep properties. However, their respective contributions on the plasticity of olivine have not been clarified yet. The potential role of cross slip will be addressed in a further study. In particular, we employ 2.5-dimensional (2.5D) dislocation dynamics simulations to investigate the interplay between thermally activated glide and climb motion and to study the effect of climb on olivine creep strain rates. Within this approach we adopt a two-dimensional (2D) reference system. Dislocations are approximated by parallel straight segments and their dynamics is followed in a reference plane. In the so called 2.5D dislocation dynamics, additional local rules

<footnote>\*francesca.boioli@univ-lille1.fr</footnote>

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 3.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
are implemented to reproduce the relevant three-dimensional dislocation mechanisms, as originally proposed to model fcc metal plasticity [5,20]. Finally, in order to describe the climb motion, the transport of matter through vacancy diffusion is taken into account and directly related to the climb dislocation velocity, similarly to Refs. [13,15]. To our knowledge, this is the first time that dislocation dynamics simulations have been employed to calculate the creep strain in a material characterized by high lattice friction, i.e., where dislocations in different slip systems glide with thermally activated mobility laws.

## II. METHOD

The description of the individual and collective behavior of dislocations is essential to determine the plastic flow of crystalline solids. Discrete dislocation dynamics is a well-established simulation technique aimed at reproducing such behavior at the mesoscale. This method is based on continuum elasticity theory that provides the description of the elastic strain field induced by dislocations, their reciprocal interaction, and the interaction of these line defects with an external stress field [10]. Moreover, in DD simulations it is possible to take into account the relevant atomistic processes controlling dislocation motion and reactions by including local rules [2,3]. Several 2D- and 3D-DD codes have been developed. 3D-DD simulations have the advantage of taking into account the topology and the evolution of curved dislocation lines in the three-dimensional crystal lattice. However 3D-DD simulations are computationally expensive and only a small amount of plastic strain can be simulated. In some cases a good insight of the relevant physical mechanisms governing dislocation motion can be achieved by using a 2D system, especially when accurate information regarding the fundamental atomistic mechanisms are lacking, as for the climb case. Here a so called 2.5D-DD code has been used, as proposed by Gómez-García et al. [5,21], where local rules have been implemented to overcome the artifacts introduced by the 2D description and to reproduce as closely as possible the 3D dislocation evolution. This code has been extended to include dislocation climb, similarly to Ref. [15]. In this section we briefly recall the main

features of the 2.5D-DD code, for more details see Refs. [5,21], and we discuss the modifications introduced in order to model olivine and to implement climb.

Olivine [(Mg,Fe)₂SiO₄] is an orthorhombic silicate with lattice parameter *a* = 4.76 Å, *b* = 10.21 Å, and *c* = 5.98 Å [22]. Numerous experiments have been carried out on both olivine (which usually contains approximately 10% of iron) and its iron-free endmember forsterite (Mg₂SiO₄) in order to characterize the slip geometry and activity in different experimental conditions [23–26]. The commonly observed dislocations have Burgers vector equal to the two shortest lattice vectors *a* [100] and *c* [001], while [010] dislocations have very scarcely been observed. Generally, a prevalence of [001] slip activity has been observed at low temperature (below 1000°C) and high differential stress, while, in the high temperature regime, plasticity is mostly dominated by [100] dislocations. Both type of defects exhibit a marked crystallographic orientation and several slip planes have been identified by different authors. In particular, [001] dislocations glide on (001), {110} , and (010) planes and are characterized by long screw segments [25–28], while [100] dislocations mostly glide on (010), (001) planes and are characterized by long straight edge segments [26,29]. At high temperature cross slip and climb are believed to play a role in olivine plasticity. Evidences of cross slip have been found above 1000°C [28,29], while climb is considered to be important above 1300°C [30]. Here we focus on the temperature range between 1400 and 1700 K, relevant for Earth’s upper mantle deformation, and on the [100](001) slip system, which is expected to dominate in the high temperature regime. In order to validate our methodology by comparing the 2.5D- with the 3D-DD simulations performed in multislip condition at low temperature [9], we performed DD simulations by including the description of a second slip system, the [001](100), expected in the low temperature regime. The reference plane employed in our simulations is (010) and the simulation area is a square of size *L*<sub>x</sub> = *L*<sub>y</sub>. In Fig. 1(a) the simulation box and the two slip systems are sketched. For the sake of simplicity we considered two sets of infinite, straight edge dislocations, perpendicular to the reference plane. Their

<figure><img src="image_2.png" /><figcaption>FIG. 1. (Color online) (a) Sketch of the 2D simulation box. Double slip conditions are used; the loading axis is parallel to the [101] direction, so that the two slip systems are homogeneously loaded. (b) Stress vs strain curves obtained by 2.5D-DD simulations performed in constant strain rate conditions. The flow stress decreases with the increasing temperature and with the decreasing strain rate values.</figcaption></figure>
<table><caption>TABLE I. Parameters used to define the glide mobility law for [100] and [001] dislocations as taken from Ref. [9].</caption><thead><tr><th></th><th>b = [100]</th><th>b = [001]</th></tr></thead><tbody><tr><td>τ<sub>P</sub> (MPa)</td><td>7170 ± 894</td><td>1730 ± 41</td></tr><tr><td>p</td><td>0.2936</td><td>0.5988</td></tr><tr><td>q</td><td>0.8475</td><td>1.1506</td></tr><tr><td>ν<sub>0</sub> (m/s)</td><td>1.65 × 10<sup>3</sup></td><td>2.75 × 10<sup>3</sup></td></tr><tr><td>l<sub>c</sub> (nm)</td><td>1.22</td><td>1.95</td></tr><tr><td>ΔH<sub>0</sub> (eV)</td><td>5.7</td><td>5.4</td></tr></tbody></table>

Burgers vector lies in the reference plane and defines the slip direction. Here we took the two slip directions inclined by 45° with respect to the x and y axis. The climb direction also lies in the reference plane and is orthogonal to the slip direction. Elastic interactions are calculated by using isotropic elasticity theory. All the simulations were performed by using the elastic constants of olivine: Poisson ratio $v = 0.25$ and shear modulus $\mu = 80$ GPa. Influence of lattice friction and interactions with point defects on dislocation velocities are taken into account via the mobility laws defined for glide and climb. In olivine, dislocations exhibit marked crystallographic orientation; this reflects a high glide resistance, probably due to an extended core structure [31]. In high lattice resistance materials, it is commonly assumed that glide is controlled by a thermally activated kink-pair mechanism. Then, the velocity of a dislocation of length $L$ can be given by an Arrhenius rate equation [32]:

$$v_g = b \frac{L}{l_c} v_D \frac{b}{l_c} \exp \left( -\frac{\Delta H(\tau^*)}{k_B T} \right), \quad (1)$$

where $b$ is used in the following to indicate the magnitude of the Burgers vector, $l_c$ is the critical length for kink nucleation, $v_D$ is the Debye frequency, $k_B$ is the Boltzmann constant, $T$ is the temperature, and $\Delta H(\tau^*)$ is the activation enthalpy of kink-pair formation which depends on the effective resolved shear stress $\tau^*$. At each dislocation position $\tau^* = \tau_{app} + \tau_{int}$ is calculated as the sum of the applied stress and the elastic interaction stress induced by all the other dislocations, projected along the slip direction. The activation enthalpy is parametrized by following the formalism of Kocks et al. [32]:

$$\Delta H(\tau^*) = \Delta H_0 \left[ 1 - \left( \frac{\tau^*}{\tau_P} \right)^p \right]^q, \quad (2)$$

where $\tau_P$ is the Peierls stress, $p$ and $q$ are empirical parameters, and $\Delta H_0$ is the total activation enthalpy. The glide mobility has been defined by using the parameters reported in Ref. [9] and shown in Table I, where the velocity for [100] and [001] dislocations has been obtained by fitting experimental data from various sources, see Ref. [9] and references therein. In 3D-DD simulations, where dislocation lines are discretized into finite segments, $L$ is effectively the length of dislocation line segments parallel to directions with high lattice friction. In a 2D framework, all dislocation lines are normal to the reference plane and supposed straight with some ideal length $L$. Here $L$ was set equal to $1 \mu$m, which is approximately the average segment length estimated from 3D-DD simulations of olivine in plastic deformation conditions closed to the one we

investigate in the present study [9]. Contrary to fcc metals, junction formation in olivine is a rare event, as predicted by Durinck et al. [9]. For this reason we eliminated the possibility for dislocations to react and form junctions. The strongest interaction expected in olivine is dipole formation. In our simulations we allowed dislocations involved in a dipole to mutually annihilate when the dipole height is smaller than a critical distance $r_a = 10b$. Our results are however rather insensitive to the value of this parameter. Hence $r_a$ can be eventually increased to speed up the simulations when a large simulation area is used. A multiplication rule is used to reproduce the general observation that the dislocation density $\rho$ increases linearly with the plastic strain $\epsilon$: $d\rho/d\epsilon = m$, where $m$ is a constant [5]. Here the value $m = 2 \times 10^{15} \text{ m}^{-2}$ has been taken to reproduce the dislocation density evolution of 3D-DD simulations [9,33]. In Fig. 1(b) we show the stress-strain curves obtained by our 2.5D-DD simulations at constant strain rate in double slip conditions. The box has been loaded along the y axis so as to induce the same Schmid's factors on the two slip systems. This corresponds to the so-called [101]$_c$ orientation in olivine where the loading axis is at 45° between [100] and [001]. The linear box size $L_x = L_y$ is $6.4 \mu\text{m}$ and the initial dislocation density is $2.6 \times 10^{12} \text{ m}^{-2}$. As we can see in Fig. 1, the flow stress increases with decreasing temperature and with increasing values of the applied strain rate. These results well reproduce the flow stress behavior obtained with 3D-DD simulations as reported Ref. [9]. The strongest interactions observed in our 2D simulations are the repulsive and attractive interactions between dislocations with the same and opposite Burgers vector, respectively. This is in agreement with the 3D case, where elastic interactions between [100] and [001] dislocations were shown to be weak [9]. Due to this result and since at elevated temperature olivine plasticity is dominated by [100] dislocations, we decided to adopt single slip conditions to address the influence of climb on high temperature plasticity.

The climb motion occurs through the absorption/emission of point defects by the dislocations. Even though both vacancies and interstitials can contribute to climb, the vacancy concentration is generally larger than the interstitial one due to their lower formation enthalpy. Here we assume that climb is controlled by vacancy diffusion and that the dislocation line is saturated with jogs. Assuming steady state conditions, the net flux of vacancy from and to the dislocation core is calculated by solving the diffusion equation and an analytical expression of the climb velocity $v_c$ can be derived [10,11]:

$$v_c = \eta \frac{D^{\text{sd}}}{b} \left[ \exp\left(\frac{\tau_c^* \Omega}{k_B T}\right) - \frac{c_{\infty}}{c_0} \right], \quad (3)$$

where $D^{\text{sd}}$ is the vacancy self-diffusion coefficient, $\Omega$ is the vacancy formation volume, and $\eta$ is a geometrical factor which depends on the geometry of the flux field. Assuming a cylindrical symmetry around the dislocation $\eta$ is given by $\frac{2\pi}{\ln(R/r_c)}$, where $R$ and $r_c$ are the radius of the two cylindrical surfaces around the dislocation core through which the vacancy flux is calculated. $r_c$ is the core radius and $R$ can be taken as half of the average dislocation distance. Since $R$ and $r_c$ appear in the logarithmic term in Eq. (3), setting $10 < R/r_c < 1000$ does not affect significantly the value of the climb velocity. Hence
we imposed $R/r_c$ to be constant and equal to 100. $c(R) = c_\infty$ is the vacancy concentration at a distance $R$ far from the dislocation, while $c_0 = \exp(-\frac{\Delta H_f}{k_B T})$ is the equilibrium vacancy concentration in a bulk at a given temperature $T$, with $\Delta H_f$ being the vacancy formation enthalpy. The self-diffusion coefficient can be also written as $D^{sd} = c_0 \exp(-\frac{\Delta H_m}{k_B T})$, with $\Delta H_m$ being the vacancy migration enthalpy. The two terms in the parentheses in Eq. (3) represent the two driving forces for climb. The first is the Peach-Kohler force in the climb direction. In fact, in the exponential term, $\tau_c^*$ is the effective stress for climb, which is given by the sum of the applied and the internal stress projected along the climb direction. Here the stress is taken to be positive when it tends to favor vacancy emission and negative when it tends to favor vacancy absorption. The second term is called a “chemical force” because it arises from a gradient in the vacancy concentration. In supersaturation conditions when $c_\infty > c_0$ this term becomes important and dislocations can climb even in absence of a mechanical force. The comparison of Eq. (3) with atomistic simulations on iron [34,35] and with phase field simulations on aluminum [36] lead to a remarkable agreement, further validating the use of such an expression to model climb with DD simulations. Unlike simple metals, vacancy diffusion in a multicomponent material involves different ionic species. Since in olivine the Si diffusion coefficient is much smaller than all the other diffusing species, we assumed that the climb process is controlled by Si diffusion. This is further justified by experiments on both pure forsterite and olivine where the measured activation enthalpy for creep and the Si self-diffusion enthalpy have comparable values [37–39]. As we can see in Eq. (3), the climb velocity depends on the average vacancy concentration in the bulk $c_\infty$. In principle, this quantity is a function of the spatial coordinates and depends on the distribution of dislocations that act as sources or sinks of vacancies. In this work we assumed that far from the dislocations the vacancy concentration is constant and equal to the equilibrium concentration in a bulk ($c_\infty = c_0$). By inserting $c_\infty = c_0$ and $D^{sd} = D_{Si}^{sd}$ in Eq. (3), we can write the climb velocity for olivine as

$$v_c = \eta \frac{D_{Si}^{sd}}{b} \left[ \exp\left(\frac{\tau_c^* \Omega}{k_B T}\right) - 1 \right], \quad (4)$$

with $\Omega = 72.4 \text{ Å}^3$, vacancy formation volume for a Shottky defect (Mg,Fe)₂SiO₄ estimated from the unit cell volume of olivine, and $D_{Si}^{sd} = D_0 \exp(-\frac{\Delta H_{Si}^{sd}}{k_B T})$, self-diffusion coefficient of silicon taken from diffusion experiments carried out on dry forsterite at ambient pressure, where the values $D_0 = 2.51 \times 10^{-7} \text{ m}^2/\text{s}$ and $\Delta H_{Si}^{sd} = 4.25 \text{ eV}$ have been measured [38]. $\Delta H_{Si}^{sd} = \Delta H_{Si}^f + \Delta H_{Si}^m$ is the activation enthalpy for Si self-diffusion which is given by the sum of the enthalpy for the formation $\Delta H_{Si}^f$ and the migration $\Delta H_{Si}^m$ of a Si vacancy. We adopted this value for the Si diffusion coefficient because it is comparable with the one measured in iron-bearing olivine, see Ref. [40], where $\Delta H_{Si}^{sd} = 4.03 \pm 0.31 \text{ eV}$ was measured, and it falls in the middle of the range of experimental values of activation enthalpy reported in literature; for example, $\Delta H_{Si}^{sd} = 5.48 \pm 0.42 \text{ eV}$ was found in Ref. [37] and $\Delta H_{Si}^{sd} = 3.02 \pm 0.16 \text{ eV}$ in Ref. [41].

<figure><img src="image_3.png" /><figcaption>FIG. 2. (Color online) Ratio of the glide over the climb velocity in olivine as a function of the temperature and the applied stress. The glide velocity is always significantly larger than the climb velocity.</figcaption></figure>

The difficulty of running simulations that simultaneously treat the glide and the climb motion arises from the different time scales that characterize the two mechanisms. In fact, by plotting the ratio between the glide and the climb velocities, $v_g$ and $v_c$, according to Eqs. (1) and (4), respectively, we can observe in Fig. 2 that glide is always two or three orders of magnitude faster than climb. It is interesting to notice that while the ratio between the glide and the climb velocity in olivine significantly increases with the increasing stress, it mildly varies with temperature. Also, the two mobility laws approach each other at very low stress and intermediate temperature, below 1400 K, as we can see in Fig. 2. This indicates that climb could play an important role even at intermediate temperatures.

In order to perform creep simulations and to resolve both glide- and climb-related events, we adopted a scheme close to the one used in Ref. [15]. First, the creep stress is applied and dislocations move by glide only using a small time step $dt_g$. When the plastic strain saturates, i.e., the dislocations reach a quasiequilibrium configuration, a larger time step $dt_c$ is used and dislocation are allowed to move by climb. In particular, the plastic strain $\Delta\epsilon$ produced at each glide step is monitored; when $\Delta\epsilon$ is lower than a critical value $\Delta\epsilon_{crit}$ over a large enough number of steps $N_{test}$ we switch to $dt_c$ and dislocations are moved by climb according to Eq. (4). When a climb displacement of one Burgers vector is achieved by at least one dislocation, the time step is switched back to $d_g$ and dislocations are relaxed by glide only. This procedure is repeated iteratively during the simulation. Here $dt_c = 10^3 dt_g$, $\Delta\epsilon_{crit} = 10^{-18}$, and $N_{test} = 400$ have been used. Periodic boundary conditions along the x and y directions have been applied in all the simulations to mimic bulk conditions and the long range stress contribution is calculated with the fast multipole algorithm. Also, we tested different ratio between the linear box sizes $L_x$ and $L_y$ and we checked that the results are independent of this ratio. Hence, the simplest solution $L_x/L_y = 1$ is used in the following section.

### III. RESULTS

The results of high temperature DD creep simulations obtained by activating the sole glide mechanism are shown
<figure><img src="image_4.png" /><figcaption>FIG. 3. (Color online) High temperature (1600 K) creep (a) and constant strain rate (b) simulations obtained by activating the glide mechanism only. (a) Plastic strain vs time curves obtained by applying a creep stress ranging from 10 to 100 MPa. (b) Stress $\sigma$ vs strain $\varepsilon$ curves obtained by applying a constant strain rate.</figcaption></figure>

in Fig. 3. All the results in this section have been obtained by imposing uniaxial loading and single slip, the [100](001), condition. The linear average dimension of the simulation box is 6.4 µm and the temperature is 1600 K. In Fig. 3(a) the strain versus time curves for different values of the creep stress, ranging from 10 to 100 MPa, are plotted. In these conditions, the plastic strain initially increases with time and then it reaches a constant value. This indicates that, after a small transient where dislocations move under the applied stress, they are blocked in a quasiequilibrium state where the glide force due to the applied stress cannot overcome the one induced by the internal stress. This interpretation is further supported by constant strain rate DD simulation results. In Fig. 3(b) the stress-strain curves obtained for different strain rate values, ranging from $10^{-4}$ to $5 \times 10^{-7}$ s$^{-1}$, are shown. For large values of the applied strain rate ($10^{-4}$, $5 \times 10^{-5}$ s$^{-1}$), after the flow stress is reached the stress remains almost constant. On the contrary, at lower strain rates (below $10^{-5}$ s$^{-1}$), an hardening behavior is observed. In the latter case, a small and almost constant strain hardening rate is found ($\sim 0.7 \times 10^{-3}$ µ). The results shown in Fig. 3 suggest that the deformation behavior, at high temperature and low applied stress (or small strain rate),

is controlled by dipolar short-range interactions that explain the formation of dislocation organized microstructures. This feature is equivalent to the forest interactions controlling dislocation microstructure formation in many other materials plasticity. Furthermore, by considering the glide mechanism only we could not reproduce the steady state creep behavior observed in experiments.

The contribution of the climb mechanism to the creep behavior of olivine is illustrated in Fig. 4, where DD simulation results for a creep stress of 40 MPa at 1600 K are presented. In particular, we compare the results obtained by activating glide and climb simultaneously with the results obtained by allowing dislocations to move by glide only. A sketch of the simulation box and of the loading condition is shown in Fig. 4(a). The tensile loading axis is inclined by 60° with respect to the glide displacement direction, which in our case coincide with the Burgers vector direction. In Fig. 4(b) the initial dislocation microstructure is shown. The initial dislocation density is $1.4 \times 10^{12}$ m$^{-2}$ and the box size is 8.3 µm. Dislocations with positive (negative) Burgers vector are represented by red (blue) full circles. Figures 4(c) and 4(d) show the final microstructure obtained when only glide is considered [Fig. 4(c)] and when both glide and climb are included in the simulation [Fig. 4(d)]. We can observe in both simulations the formation of dislocation walls perpendicular to the glide direction, as dislocations with the same sign tend to self-organize along this direction. Moreover, most of the dislocations are trapped into dipole configurations, as we can observe in Figs. 4(c) and 4(d), where we see numerous couples of dislocations with opposite signs closely packed together. In Fig. 4(e) the strain vs time curve obtained by including the climb mechanism (blue curve) is compared to the curve obtained by activating the glide mechanism only (red curve). When climb is active a steady state regime is found and the plastic strain increases linearly with time. In Fig. 4(f) the contribution of glide (green) and of climb (red) to the total plastic strain (blue) is shown. As we can see, even if climb is active, the plastic strain is almost completely produced by dislocations moving in their glide planes. The effect of climb is to release dislocations from a “jammed” configuration, i.e., when dislocations are trapped in a minimum energy configuration, so that they can further move by glide and produce plastic strain. At a given time most of the dislocations are immobile and stored in a dislocation microstructure. Climb allows dislocations involved into dipoles to move outside the glide plane to annihilate or to bypass repulsive dislocations. As a consequence, the dislocation microstructure is relaxed and other dislocations become free to move by glide until a new jammed configuration is reached. The repetition of this process in time leads to the steady state creep condition, characterized by a linear increase of the plastic strain with time and a nearly constant dislocation density. In Fig. 5 DD creep curves obtained by varying the applied stress from 10 to 60 MPa at 1600 K are shown. In all simulations we observe a linear dependence of the plastic strain on time and larger strain rates are found for larger values of the creep stress, as can be seen in Fig. 5(a). Moreover, after an initial transient where dislocation multiplication is observed, the dislocation density fluctuates around a constant value, see Fig. 5(b). In this stage both dislocation annihilation and multiplication operate
<figure><img src="image_5.png" /><figcaption>FIG. 4. (Color online) (a) Sketch of the 2D simulation box. Single slip conditions are used. The angle between the loading axis and glide displacement direction is 60°. (b) Initial dislocation microstructure and final dislocation microstructure in glide only (c) and glide coupled with climb (d) dislocation motion. (e) Plastic deformation ε vs time, as obtained in creep conditions, assuming glide only (red) and glide coupled with climb (blue) dislocation motion. (f) The plastic strain produced by glide ε$_{glide}$ (green curve) and climb motion ε$_{climb}$ (red curve) as a function of time is compared with the total plastic strain ε$_{tot}$ (blue curve).</figcaption></figure>

so as to maintain an equilibrium dislocation density value. The latter increases with the increasing creep stress values as dislocation dipole height can be smaller at higher stress. From our DD simulations it emerges that deformation involving glide assisted by climb leads to a steady state creep behavior. Moreover, in order to describe at least qualitatively plastic deformation of olivine at high temperature, it is necessary to take climb into account.

A large number of DD simulations, taking into account both the glide and climb, have been performed in order to calculate the high temperature creep strain rate in olivine. The range of applied stresses σ and temperatures T are, respectively, from 10 to 115 MPa and from 1400 to 1700 K. In particular, for each couple of values (σ,T) we run numerous simulations varying the average linear box size, from 3.8 to 21.8 µm, and the initial dislocation configuration, so as to avoid any influence of the boundary condition and of the initial dislocation microstructure on the calculated creep rates. The latter were determined by a linear least-square fit of the plastic strain vs time curves and then by taking the average value for each set of simulations run at constant σ and T. Similarly, we calculated the average dislocation density. In Fig. 6 we plot the average dislocation density as a function of the shear stress resolved in the glide plane τ. We observe that the average dislocation density does not depend

on temperature. This result confirms that the flow stress is controlled by short range dislocation-dislocation interactions in the dislocation microstructure. In view of the discussion in Sec. IV about the mobile dislocation density and velocity, we analyze the dependence of the total dislocation density on the resolved shear stress τ. Typically, when dislocation-dislocation interactions govern the microstructure and the flow stress, it is observed that the dislocation density is proportional to the square of the critical flow stress and can be described by rewriting Taylor’s equation as follows:

$$ \rho = \left( \frac{\tau}{\mu} \right)^2 C/b^2. \qquad (5) $$

In Taylor’s equation the constant C is given by $1/α^2$ where α is a dimensionless parameter that reflects the strength of the dislocation-dislocation interaction, i.e., the strength of the obstacles. By fitting our DD results with Eq. (5) (solid gray curve in Fig. 6) we find $C = 21.6 \pm 0.7$ and consequently $α = 0.21$. This value of α is consistent with the case where the microstructure is controlled by dipole interactions [42]. However, we can see in Fig. 6 that Eq. (5) does not reproduce very well the evolution of the dislocation density at low stresses. Indeed Taylor’s equation usually holds for materials with low lattice resistance when the main
<figure><img src="image_6.png" /><figcaption>FIG. 5. (Color online) (a) Creep curves of the total plastic strain vs time for different values of the applied stress at T = 1600 K. (b) Corresponding evolution of the dislocation density with time.</figcaption></figure>

strengthening mechanism is forest interactions; for example it applies remarkably well for fcc metals in a wide range of applied stresses [43]. In olivine, at low and moderate temperature, large α coefficient >1, which are not consistent

<figure><img src="image_7.png" /><figcaption>FIG. 6. (Color online) Total dislocation density ρ as a function of the applied stress resolved in the glide plane τ.</figcaption></figure>

with a dislocation interaction strengthening mechanism, are generally found [9,23]. This can be explained by the high lattice friction of this material, for which dislocations cannot easily bow out to bypass obstacles. At high temperature, climb assists dislocation motion and helps them bypassing obstacles. Still, dislocations exhibit marked crystallographic orientation which reveal that the contribution of lattice resistance to the material strengthening is not negligible even in the high temperature regime. The discrepancy between the simulation results and Eq. (5) at low stress also points in this direction. If we analyze the relationship between the dislocation density and the applied stress with the more general equation:

$$ \rho = \left(\frac{\tau}{\mu}\right)^s C'/b^2, \qquad (6) $$

we find that the best fit for the stress exponent s is 1.31 ± 0.03, with C' = 0.104 ± 0.0025. This curve, shown with a dotted black line in Fig. 6, reproduces well the simulation results. Moreover, it is in good agreement with several experimental findings. In Ref. [29] a power law with a stress exponent s = 1.37 has been found by fitting the values of the dislocation density, measured in olivine samples deformed at high temperature and low strain rates, as a function of the stress. Similarly, s = 1.41 ± 0.16 was found by analyzing olivine samples deformed at high temperature and pressure of 1–2 GPa [44].

High temperature deformation is commonly analyzed by using the following equation for the creep rates:

$$ \dot{\epsilon} = \dot{\epsilon}_0 \sigma^n \exp \left( -\frac{Q}{k_B T} \right), \qquad (7) $$

where Q is the creep activation enthalpy, σ is the applied stress, and T is the temperature. Q, n, and $\dot{\epsilon}_0$ are generally assumed to be constant and depend on the physical mechanism that control the plastic behavior. In Fig. 7(a) the values of the strain rates as obtained by DD simulations are reported as a function of the reciprocal temperature T. From the negative slope of the logarithm of the strain rate versus the reciprocal temperature curves, it is possible to extract the activation enthalpy Q. Experimentally this value has been found to be close to the Si self-diffusion activation enthalpy $\Delta H_{Si}^{sd}$. In agreement with the experiments [23,24,26,45], we find Q = 5.07 eV, value which is in between the activation enthalpy for Si self-diffusion ($\Delta H_{Si}^{sd} = 4.25$ eV) and kink-pair activation enthalpy for [100] dislocations ($\Delta H_0 = 5.7$ eV). This indicates that dislocation mobility in olivine is controlled by the interplay between two kinetics, the one of glide and the one of climb. The stress exponent n in Eq. (7) can be calculated as the derivative of the logarithm of the strain rate with respect to the logarithm of stress. In Fig. 7(b) the creep strain rates are reported as a function of the applied stress for four different T values. A constant slope of the log $\dot{\epsilon}$ curve is found for all temperatures considered. This corresponds to a constant stress exponent value n close to 3 in the (σ, T) range here explored. In particular n is equal to 3.1, 3.2, 3.2, and 2.9 at 1700, 1600, 1500, and 1400 K, respectively. In Sec. IV we compare our results with olivine relevant creep models and we propose a semianalytical model based on our DD simulation results aimed at describing
<figure><img src="image_8.png" /><figcaption>FIG. 7. (Color online) DD strain rate values $\epsilon'$ as a function (a) of the reciprocal temperature $T$ and (b) of the applied creep stress $\sigma$.</figcaption></figure>

the creep strain rates in olivine and possibly in other high lattice friction materials.

## IV. COMPARISON WITH CREEP MODELS AND DISCUSSION

Several efforts have been made in the past to define simple analytical equations that can describe the creep strain rates in different deformation conditions. Most of the theoretical works on dislocation creep are focused on fcc metals and are based on the pioneer works of Weertman [46,47] and Nabarro [48]. These works helped explain the steady state creep behavior in some limiting cases, but a comprehensive creep model is still lacking. First, it is generally assumed that the creep strain rates can be written in terms of Orowan's equation:

$$ \dot{\epsilon} = \rho_m b v_m, \quad (8) $$

where $\rho_m$ and $v_m$ are the average mobile dislocation density and velocity. Then, few assumptions are made to evaluate the latter two quantities $\rho_m$ and $v_m$. In most of the cases, $\rho_m$ is substituted with the total dislocation density $\rho$ and $v_m$ with the instantaneous glide or climb dislocation velocity $v$. In climb-controlled models, it is assumed that climb is much slower than glide. Hence, the time $t_g$ needed to glide a distance $L$ is much smaller than the time needed to climb from one

glide plane to another one. The dislocation velocity, written as $v = \frac{L}{t_g+t_c}$, is then simplified as $v = \frac{L}{t_g+t_c} \sim \frac{L}{t_c} = \frac{L}{d} v_c$, where $d$ is the distance between glide planes and $v_c$ is the climb velocity. By substituting this expression in Eq. (8), it results in $\dot{\epsilon} = \frac{L}{d} \rho b v_c$. By further assuming that the dislocation density is proportional to the square of the applied stress, according to the Taylor relationship [Eq. (5)], and that, for small stress values, the climb velocity varies linearly with stress [when $\tau \Omega \ll k_B T$ it is possible to simplify $\exp(-\frac{\tau \Omega}{k_B T}) - 1 \sim \frac{\tau \Omega}{k_B T}$ in Eq. (3)], the strain rates can be described by a power law with a stress exponent equal to 3. A creep power law characterized by the same stress exponent has been found in the model of Nabarro [48], where creep is produced by dislocation climb only. Similarly to the climb-controlled model it is assumed $\rho \propto \tau^2$ and that the climb velocity varies linearly with the stress, for small stress values. The dislocation velocity is then assumed to be exactly the climb velocity, obtaining $\dot{\epsilon} \propto \tau^3$. Stress exponent larger than 3, as found for fcc Al, could be reproduced by making some ad hoc assumptions on the multiplication and annihilation processes [47,49]. Only recently, DD simulations on Al micropillars have been used to describe the different stress exponents, ranging from 1 to 6, in different stress and temperature conditions by coupling the glide and the climb mechanisms [15]. This shows the capabilities of DD simulations to predict the collective behavior of dislocations. Still, in olivine, a power law characterized by a stress exponent close to 3–3.5 is found in a wide range of temperature and stress conditions [50,51]. Our DD results are in good agreement with these experimental evidences. In fact, we were able to capture the variation of both the dislocation density and the creep strain rate with the applied stress and temperature. In the following we compare our numerical results with the theoretical models proposed in literature and we derive a semianalytical model based on the analysis of our dislocation dynamics simulations.

In order to make use of Eq. (8) to extract a simple expression for the creep strain rates, the first issue is to verify that the collective behavior of dislocations can be described as an average of the single dislocation behavior. Then we need to evaluate the mobile dislocation density and velocity that characterize the collective dislocation motion. To this aim we analyze our DD simulation results and extract the average mobile dislocation density $\rho_m$ and velocity $v_m$ for the ensemble of dislocations that contribute, respectively, to 99% and 95% of the total plastic strain produced during the simulation $\epsilon_{tot}$. Details about how we calculate the mobile dislocation velocity $v_{m,99}$ ($v_{m,95}$) and density $\rho_{m,99}$ ($\rho_{m,95}$) are reported in the Appendix. The results of the analysis of our DD simulations performed at 1400, 1600, and 1700 K, are shown in Figs. 8 and 9. The average density of mobile dislocations obtained for the two criteria are plotted in Fig. 8(a) and compared with the total dislocation density $\rho_{tot}$. Most dislocations are immobile. Figure 8(b) shows that only a small fraction of dislocations effectively contribute to the plastic strain. Almost all the plastic strain (99%) produced during creep results from the motion of a small portion of the dislocations, between 30% and 40%, approximately. Among those, a smaller ensemble of mobile dislocations, 10% of the total, is responsible for 95% of the total plastic deformation. On one hand, this fraction is almost constant and does not significantly change with the applied stress or temperature. Interestingly we found
<figure><img src="image_9.png" /><figcaption>FIG. 8. (Color online) (a) Dislocation density as a function of the resolved stress $\tau$. Solid (dashed) lines are used to plot the total (mobile) dislocation density. (b) Fraction of mobile dislocation as a function of the resolved stress $\tau$. In both panels red diamonds, green boxes, and blue circles are used to display data at 1700, 1600, and 1400 K, respectively. Open (filled) symbols are used to indicate the average mobile dislocation density that contributes to 99% (95%) of the total plastic strain $\epsilon_{tot}$.</figcaption></figure>

that also in constant strain rate condition plastic deformation is sustained by approximately the same constant fraction of mobile dislocations. On the other hand, the velocity of mobile dislocations strongly depends on temperature, since both glide and climb processes are thermally activated. In Figs. 9(a)–9(c) the values obtained for the mobile dislocation velocity $v_{m,99}$ ($v_{m,95}$) are plotted as a function of the applied stress and compared with the instantaneous glide and climb velocities. Figure 9 demonstrates that the mobile dislocation velocity cannot be replaced by either the pure climb or the climb velocity. It clearly results from an interplay between the two kinetics. In order to derive an expression for the mobile dislocation velocity, we can make few assumptions based on our DD results. First, we can assume that the distance covered by climb $\Delta l_c$ is negligible compared to the glide mean free path $\Delta l_g$ traveled by mobile dislocations and express the distance $\Delta l$ traveled in the time interval $\Delta t$, as the glide mean free path $\Delta l \cong \Delta l_g = v_g(\tau, T)\Delta t_g$. This is supported

by the results in Fig. 4(f), where it emerges that the entire plastic strain is only produced by glide during the whole creep simulation. On the contrary, the waiting time $\Delta t_c$ needed to release a dislocation from a quasiequilibrium configuration is much larger than the time $\Delta t_g$ needed to glide from a jammed configuration to another one, so that $\Delta t = \Delta t_c + \Delta t_g \cong \Delta t_c$. This is supported first by the observation that the glide velocity is two or three order of magnitude larger than the climb velocity and it is corroborated by the observation that a large fraction of dislocations is nearly immobile. By following these two considerations, we can write the average mobile dislocation velocity as

$$v_m = \frac{\Delta l}{\Delta t} \cong \frac{v_g(\tau, T)\Delta t_g}{\Delta t_c + \Delta t_g} \cong v_g(\tau, T) \frac{\Delta t_g}{\Delta t_c}. \quad (9)$$

Here the mobile dislocation velocity is expressed as the product of the pure glide velocity and the ratio between the time needed to glide the mean free path and the waiting time between two glide displacements. We notice that $\Delta t_c$ is not simply the time needed to do a climb step, but it represents the time interval during which the dislocation is trapped in a local equilibrium configuration. $\Delta t_c$ depends on the local dislocation configuration: a dislocation remains immobile until it climbs to bypass an obstacle or until one or few of the neighboring dislocations move or annihilate by climb so as to modify the local internal stress field and release it. In general, the instantaneous glide velocity for each dislocation is due to an effective stress $\tau^*$ given by the sum of the applied stress $\tau$ and of the internal stress $\tau_{int}$. While for the immobile dislocations $\tau \sim -\tau_{int}$, for the subset of mobile dislocations we can assume that the average stress value controlling the glide process is given by the applied creep stress $\tau > \tau_{int}$. In Fig. 9(d) we plot the ratio between the average mobile dislocation velocity $v_m$ and the glide velocity $v_g$ as a function of $\tau$ at three different temperatures $T$. This ratio turns out to be approximately constant with the increasing stress values. By setting $k = v_m(\tau, T)/v_g(\tau, T)$, we fit the ratio between the mobile dislocation velocity $v_{m,99}$ ($v_{m,95}$) and the glide velocity and we extract the $k$ value for each curve of $v_m$. The curves obtained by this fitting are plotted in Figs. 9(a)–9(c), dashed lines in the graphs. The values of $k$ extracted from the fit of $v_{m,99}$ and $v_{m,95}$ are, respectively, 0.019 and 0.074 at 1700 K, 0.021 and 0.078 at 1600 K, and 0.010 and 0.034 at 1400 K. Since the values of $v_{m,99}$ and $v_{m,95}$ calculated from our DD simulations are well described by the fitting curves, we can infer that the mobile dislocation velocity is given by the glide velocity reduced by a constant $k$. By using Eq. (9) the latter can also be interpreted as the ratio between the time needed for glide and the waiting time (controlled by climb) necessary to bypass obstacles. Hence we observe that this ratio is rather temperature insensitive. Values obtained at 1600 and 1700 K are nearly identical, while at 1400 K they are reduced by a factor of 2 only. Since $k$ represents the competition between glide and climb kinetics, we can qualitatively explain such a result from Fig. 2 which shows that the ratio between glide and climb velocities does not vary significantly with temperature (contour lines in Fig. 2 are almost vertical). In fact, if we select a given stress value in Fig. 2, we can observe that the difference between the ratio $v_g/v_c$ at 1700 and at 1400 K is
<figure><img src="image_10.png" /></figure>

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

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

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

FIG. 9. (Color online) Dislocation velocity as a function of the resolved stress $\tau$ at 1400 K (a), 1600 K (b), and 1700 K (c). The solid (dotted) curves represent the pure glide (climb) velocity. The average velocity for the subset of dislocations that contribute to 99% (95%) of the total plastic strain $\epsilon_{tot}$ are plotted with filled (open) symbols. (d) Ratio between the average mobile dislocation velocity and the glide velocity. Circle, box, and diamond symbols are used to show the results at 1400, 1600, and 1700 K, respectively. A constant ratio between the glide velocity and the mobile dislocation velocity $v_m$ is found.

only a factor of 2–3. In the examined stress range, a power law with a constant stress exponent $m$ close to 2 is found to fit the glide velocity and, consequently, the mobile dislocation velocity ($v_m \propto \tau^m$). In particular, the best fit of $m$ is 1.88, 1.85, and 1.5 at 1700, 1600, and 1400 K, respectively. It is worth to say that, while in the limited stress and temperature range considered here it is appropriate to describe the glide mobility law with a power law characterized by a constant exponent $m$, this approximation may not be valid in general.

Once having calculated the average mobile dislocation density and velocity, we can derive the strain rates from Orowan's equation [Eq. (8)]. In Fig. 10 we compare the strain rates directly obtained from our DD simulations with the ones calculated by inserting the values of $\rho_{m,99}$ and $v_{m,99}$ or $\rho_{m,95}$ and $v_{m,95}$ in Orowan's equation. The agreement is quite remarkable. This justifies the use of Orowan's equation, as long as the average mobile dislocation density and velocity are adopted. This demonstrates that the collective behavior of dislocations in creep conditions can be well described by the motion of small set of mobile dislocations. These dislocations move essentially by glide and their average velocity is reduced by a factor that depends on the time needed to climb and bypass

obstacles. In the examined range, a creep power law with a constant stress exponent close to 3 is found as the mobile dislocation density is a constant fraction of the total density which scales to $\tau^s$, with $s \sim 1.3$, and the mobile dislocation velocity increases with the creep stress as $\tau^m$, with $m \sim 1.7$.

V. CONCLUSIONS

Dislocation dynamics simulations performed with our 2.5D approach are shown to be a useful tool to model the high temperature deformation behavior in olivine. From our results it emerges that it is fundamental to consider the climb mechanism in order to reach steady state creep conditions. In the temperature range between 1400 and 1700 K and with applied creep stresses between 10 and 100 MPa, it is possible to describe the creep strain rates by a power law. A constant value for the stress exponent $n$ close to 3 and an activation enthalpy of 5.07 eV are found in agreement with published experimental results. From an analysis of the average mobile dislocation density and velocity we formulate a semianalytical model able to reproduce the creep behavior in olivine at low stress and high temperature. In particular, we
<figure><img src="image_14.png" /><figcaption>FIG. 10. (Color online) Creep strain rates as obtained by DD strain vs time curves (solid lines) and by inserting $\rho_m$ and $v_m$ in Orowan's equations. The value calculated at 1700, 1600, and 1400 K are plotted in the figure with diamond, box, and circles, respectively. The strain rates obtained by using the average mobile dislocation density and velocity for the subset of dislocations that contribute to 99% (95%) of the total plastic strain $\epsilon_{tot}$ are plotted with filled black (open gray) symbols.</figcaption></figure>

find that the collective behavior of dislocations in this regime can be described from the motion of a small ensemble of mobile dislocations that are responsible for the whole plastic strain. These dislocations are a constant fraction of the total dislocation density and have an average velocity which is given by the glide velocity reduced by a coefficient that is rather insensitive to temperature. As a consequence, even though climb is fundamental to maintain an equilibrium between the multiplication and annihilation processes, it emerges that the glide kinetics plays a major role in determining the creep strain rates and the constant exponent $n \sim 3$ in the creep power law. In olivine, even at high temperature, lattice friction is significantly large. Therefore, the glide mobility law differs substantially from the one characterizing fcc metals like Al, where a stress exponent ranging from 1 to 6 has been found by DD simulations. In the latter case, only climb is a thermally activated process and the competition between the two kinetics differs substantially at different temperatures.

Our capability to reproduce the qualitative creep behavior and the good agreement we found between the predicted and the measured stress exponent values in a wide range of temperature and stress, demonstrate the strength of our simulations. We believe that this numerical approach can be useful to investigate high temperature plasticity in other relevant minerals and to compare the creep behavior in different conditions and for different material parameters, such as the diffusion coefficient. Finally, we think that the creep model here proposed is useful to analyze the creep behavior of other materials characterized by high lattice resistance, such as bcc metals or zirconium alloys, and also to understand the influence of the glide and the climb motion on the mechanical behavior in different temperature and stress conditions.

## ACKNOWLEDGMENT

Financial support by the European Research Council under the Seventh Framework Program (FP7), ERC Grant No. 290424 Rheoman, is gratefully acknowledged.

## APPENDIX: CALCULATION OF THE AVERAGE MOBILE DISLOCATION DENSITY AND VELOCITY

In order to evaluate the mobile dislocation density $\rho_m$ and velocity $v_m$ we analyze our DD simulation results and extract $\rho_m$ and $v_m$ from the same set of simulations from which we obtained the creep strain rates and equilibrium dislocation density reported in Figs. 6 and 7. To define the ensemble of mobile dislocations we need to establish a criterion. Here we chose to select the dislocations that mostly contribute to the plastic strain $\epsilon$. In particular, we adopt two criteria: for each simulation we calculate the average dislocation density and velocity for the ensemble of dislocations that contribute, respectively, to 99% and 95% of the total plastic $\epsilon_{tot}$ strain produced during the simulation. To calculate $\rho_m$ and $v_m$ we proceed in the following way. For each simulation, first we pick a test value for the dislocation velocity $v_{test}$. For all the dislocations $N_i$ in the simulation box we calculate the distance $\Delta r_i$ traveled by the $i$th dislocation in a time interval $\Delta t \gg dt_c$ and we select the subset of $N_j$ dislocations characterized by a velocity $v_i = \frac{\Delta r_i}{\Delta t} \ge v_{test}$. For this subset we calculate the average dislocation velocity:

$$ \bar{v}_{\Delta t} = \frac{1}{N_j} \sum_{j=1}^{N_j} v_j $$

and the plastic strain produced in the time interval $\Delta t$:

$$ \Delta \epsilon_{\Delta t} = \int_{0}^{\Delta t} \sum_{j=1}^{N_j} b_j v'_j dt. $$

Here $b_j$ and $v'_j$ are the Burgers vector and the instantaneous dislocation velocity calculated at each glide or climb time step $dt$ in the interval $\Delta t$ for the $j$th dislocation. We perform this computation for subsequent time intervals $\Delta t$ until the end of the simulation ($t = t_f$) and we extract the average dislocation velocity

$$ \bar{v}_{test} = \frac{\Delta t}{t_f} \sum \bar{v}_{\Delta t} $$

and density

$$ \bar{\rho}_{test} = \frac{1}{L_x L_y} \frac{\Delta t}{t_f} \sum \bar{N}_j $$

for the subset of dislocations with $v_j \ge v_{test}$ and we calculate the total amount of plastic strained $\epsilon_{test} = \sum \bar{\epsilon}_{\Delta t}$ originated by such set of dislocations. By iterating these steps for increasing values of $v_{test}$, we defined $\rho_{m,99} = \bar{\rho}_{test}$ and $v_{m,99} = \bar{v}_{test}$ ($\rho_{m,95} = \bar{\rho}_{test}$ and $v_{m,95} = \bar{v}_{test}$) when we find $\epsilon_{test} = 0.99\epsilon_{tot}$ ($\epsilon_{test} = 0.95\epsilon_{tot}$), being $\epsilon_{tot}$ the total plastic strain generated during the simulation. For each simulation we checked different values of $\Delta t$, between $10^3 dt_g$ and $10^6 dt_g$, obtaining consistent results for the different values.
[1] P. Cordier, J. Amodeo, and P. Carrez, *Nature* (London) **481**, 177 (2012).

[2] B. Devincre, R. Madec, G. Monnet, S. Queyreau, R. Gatti, and L. Kubin, in *Mechanics of Nano-objects* (Presses de l'Ecole des Mines de Paris, Paris, 2011), p. 85.

[3] A. Arsenlis, W. Cai, M. Tang, M. Rhee, T. Oppelstrup, G. Hommes, T. G. Pierce, and V. V. Bulatov, *Model. Simul. Mater. Sci. Eng.* **15**, 553 (2007).

[4] R. Madec, B. Devincre, and L. P. Kubin, *Phys. Rev. Lett.* **89**, 255508 (2002).

[5] D. Gómez-García, B. Devincre, and L. P. Kubin, *Phys. Rev. Lett.* **96**, 125503 (2006).

[6] K. W. Schwarz, *Phys. Rev. Lett.* **91**, 145503 (2003).

[7] M. Grydlik, F. Boioli, H. Groiss, R. Gatti, M. Brehm, F. Montalenti, B. Devincre, F. Schaffler, and L. Miglio, *Appl. Phys. Lett.* **101**, 013119 (2012).

[8] J. Amodeo, P. Carrez, B. Devincre, and P. Cordier, *Acta Mater.* **59**, 2291 (2011).

[9] J. Durinck, B. Devincre, L. P. Kubin, and P. Cordier, *Am. Mineral.* **92**, 1346 (2007).

[10] J. P. Hirth and J. Lothe, *Theory of Dislocations* (Krieger, Malabar, 1992).

[11] J. L. Martin and D. Caillard, *Thermally Activated Mechanisms in Crystal Plasticity* (Pergamon, New York, 2003).

[12] N. M. Ghoniem, S.-H. Tong, and L. Z. Sun, *Phys. Rev. B* **61**, 913 (2000).

[13] D. Mordehai, E. Clouet, M. Fivel, and M. Verdier, *Philos. Mag.* **88**, 899 (2008).

[14] B. Bakó, E. Clouet, L. M. Dupuy, and M. Blétry, *Philos. Mag.* **91**, 3173 (2011).

[15] S. M. Keralavarma, T. Cagin, A. Arsenlis, and A. A. Benzerga, *Phys. Rev. Lett.* **109**, 265504 (2012).

[16] K. M. Davoudi, L. Nicola, and J. J. Vlassak, *J. Appl. Phys.* **111**, 103522 (2012).

[17] K. Danas and V. S. Deshpande, *Model. Simul. Mater. Sci. Eng.* **21**, 045008 (2013).

[18] M. Huang, Z. Li, and J. Tong, *Int. J. Plasticity* **61**, 112 (2014).

[19] G. Po and N. Ghoniem, *J. Mech. Phys. Solids* **66**, 103 (2014).

[20] A. A. Benzerga, Y. Bréchet, A. Needleman, and E. V. der Giessen, *Model. Simul. Mater. Sci. Eng.* **12**, 159 (2004).

[21] D. Gómez-García, B. Devincre, and L. P. Kubin, *MRS Proc.* **578**, 131 (1999).

[22] J. R. Smyth and R. M. Hazen, *Am. Mineral.* **58**, 588 (1973).

[23] D. L. Kohlstedt and C. Goetze, *J. Geophys. Res.* **79**, 2045 (1974).

[24] W. B. Durham and C. Goetze, *J. Geophys. Res.* **82**, 5737 (1977).

[25] M. Darot and Y. Gueguen, *J. Geophys. Res.: Solid Earth* **86**, 6219 (1981).

[26] Y. Gueguen and M. Darot, *Phys. Earth Planetary Interiors* **29**, 51 (1982).

[27] C. B. Raleigh, *J. Geophys. Res.* **73**, 5391 (1968).

[28] P. Phakey, G. Dollinger, and J. Christie, in *Flow and Fracture of Rocks* (American Geophysical Union, Washington, DC, 1972), Vol. 16, pp. 117–138.

[29] Q. Bai and D. L. Kohlstedt, *Tectonophysics* **206**, 1 (1992).

[30] B. Evans and C. Goetze, *J. Geophys. Res.: Solid Earth* **84**, 5505 (1979).

[31] P. Carrez, A. M. Walker, A. Metsue, and P. Cordier, *Philos. Mag.* **88**, 2477 (2008).

[32] U. Kocks, A. Argon, M. Ashby, B. Chalmers, J. Christian, and T. Massalski, *Progress in Material Science* (Pergamon, Oxford, UK, 1975).

[33] J. Durinck, Ph.D. thesis, University of Lille1, 2005.

[34] E. Clouet, *Phys. Rev. B* **84**, 092106 (2011).

[35] M. Kabir, T. T. Lau, D. Rodney, S. Yip, and K. J. Van Vliet, *Phys. Rev. Lett.* **105**, 095501 (2010).

[36] P.-A. Geslin, B. Appolaire, and A. Finel, *Appl. Phys. Lett.* **104**, 011903 (2014).

[37] R. Dohmen, S. Chakraborty, and H.-W. Becker, *Geophys. Res. Lett.* **29**, 24 (2002).

[38] H. Fei, C. Hegoda, D. Yamazaki, M. Wiedenbeck, H. Yurimoto, S. Shcheka, and T. Katsura, *Earth Planetary Sci. Lett.* **345-348**, 95 (2012).

[39] H. Fei, M. Wiedenbeck, D. Yamazaki, and T. Katsura, *Nature (London)* **498**, 213 (2013).

[40] S. Karato and M. Ogawa, *Phys. Earth Planetary Interiors* **28**, 102 (1982).

[41] B. Houlier, M. Cheraghmakani, and O. Jaoul, *Phys. Earth Planetary Interiors* **62**, 329 (1990).

[42] P. Haasen, in *Dislocation Dynamics in the Diamond Structure* (McGraw Hill, New York, 1968), pp. 701 and 718.

[43] M. Sauzay and L. P. Kubin, *Progress Mater. Sci.* **56**, 725 (2011).

[44] S.-I. Karato and H. Jung, *Philos. Mag.* **83**, 401 (2003).

[45] S. Mei and D. L. Kohlstedt, *J. Geophys. Res.: Solid Earth* **105**, 21471 (2000).

[46] J. Weertman, *J. Appl. Phys.* **28**, 1185 (1957).

[47] J. Weertman, *J. Appl. Phys.* **26**, 1213 (1955).

[48] F. R. N. Nabarro, *Philos. Mag.* **16**, 231 (1967).

[49] J. Weertman and J. R. Weertman, *Annu. Rev. Earth Planetary Sci.* **3**, 293 (1975).

[50] Z. M. Jin, Q. Bai, and D. L. Kohlstedt, *Phys. Earth Planetary Interiors* **82**, 55 (1994).

[51] Q. Bai, S. J. Mackwell, and D. L. Kohlstedt, *J. Geophys. Res.: Solid Earth* **96**, 2441 (1991).