A stochastic solver based on the residence time algorithm for crystal plasticity models

The deformation of crystalline materials by dislocation motion takes place in discrete amounts determined by the Burgers vector. Dislocations may move individually or in bundles, potentially giving rise to intermittent slip. This confers plastic deformation with a certain degree of variability that can be interpreted as being caused by stochastic fluctuations in dislocation behavior. However, crystal plasticity (CP) models are almost always formulated in a continuum sense, assuming that fluctuations average out over large material volumes and/or cancel out due to multi-slip contributions. Nevertheless, plastic fluctuations are known to be important in confined volumes at or below the micron scale, at high temperatures, and under low strain rate/stress deformation conditions. Here, we develop a stochastic solver for CP models based on the residence-time algorithm that naturally captures plastic fluctuations by sampling among the set of active slip systems in the crystal. The method solves the evolution equations of explicit CP formulations, which are recast as stochastic ordinary differential equations and integrated discretely in time. The stochastic CP model is numerically stable by design and naturally breaks the symmetry of plastic slip by sampling among the active plastic shear rates with the correct probability. This can lead to phenomena such as intermittent slip or plastic localization without adding external symmetry-breaking operations to the model. The method is applied to body-centered cubic tungsten single crystals under a variety of temperatures, loading orientations, and imposed strain rates.


Introduction
Plastic deformation in single crystals can generally be construed as a process in which slip on a set of well-defined crystallographic planes results in irreversible shape changes in a material [1,2]. Through its various forms, the crystal plasticity method (CP) has proven extremely successful in predicting crystalline materials deformation thanks to (i) the existence of a robust mathematical theory [1,[3][4][5], (ii) a strong connection between the physics of slip and the geometry of deformation [6], and (iii) the development of accurate and efficient integration algorithms [7][8][9][10]. While some challenges still remain, such as the indetermination of the active slip systems under a prescribed set of conditions [11], or the distinction between stored and geometrically-necessary dislocation density evolutions [12,13], suitable approximations exist that are generally valid and give accurate results in most applications.
However, the standard implementation in CP models of crystallographic slip as a linear combination of plastic shear on a finite set of slip systems assumes that these slip systems are independent of each other [14][15][16]. Such notion is generally grounded on an adiabatic interpretation of dislocation-mediated shear 1 , which is often questionable and leads to an homogeneous plastic response even in cases where valid localized deformation pathways exist [12,17]. Such homogenization is mathematically not problematic, as the underlying constitutive models intrinsically ensure stability in the solutions. However, it is often at odds physically with the observed material response [18][19][20] in which symmetry in the evolution of plastic slip can naturally break and lead to a non-homogeneous distribution of plastic shear. Accounting for symmetry-breaking modalities in the slip response has typically been done by adding external probabilistic descriptors of certain internal variables [21][22][23][24][25]. Alternatively, stochasticity can be introduced intrinsically, via the integration algorithm itself, a practice widely used for integrating differential equation systems [26][27][28][29]. As such, stochastic fluctuations are an intrinsic part of the solution, not just an externally added variable. When suitably adapted for CP models, such solvers could lead to fluctuations in plastic shear that can numerically break the symmetry in slip systems with identical Schmid factors. These fluctuations can also be regarded as reflecting the natural intrinsic variability of certain internal variables, i.e., they can be justified as 'physical' in many cases.
The main objective of this paper is to develop a stochastic solver based on the residence time algorithm [30,31] for standard single crystal plasticity problems. Under our approach, the mathematical formulation of the CP problem remains unchanged, which makes the algorithm completely general and independent of the constitutive framework and material model chosen. For proof of principle, here we focus on single-point material cases with no associated boundaryvalue problem. The paper is organized as follows. First, in Sect. 2 we provide an essential review of the theory and methods used here, including brief descriptions of the residence time algorithm and crystal plasticity approach. In Sect. 3, we present results of CP simulations of uniaxial deformation in single crystal tungsten using the proposed stochastic algorithm under a variety of different scenarios. Section 4 contains a detailed discussion of the main results as well as of the potential uses of the method. We finalize with the conclusions and acknowledgements.

The residence time algorithm
The starting point is a (continuous) transition function P i j (t) representing the probability that a random variable ξ 1 (t) 1 The process is adiabatic in the sense that, during a small time interval, slip in a given slip system is assumed to evolve without being affected by slip in the remaining slip systems, i.e., the time step is sufficiently short to justify a lack of interaction among the operative slip systems. describing the state of a continuous-time Markov chain is X j at time t having originated from an initial state X i at time t = 0, i.e. 2 : As such, all the p i (t) = j =i P i j (t) conform a vector p of probabilities of dimension N equal to the total number of states accessible to the system. A square positive-definite N × N matrix Q = {q i j } satisfying j =i q i j = r i (where r i is the exit rate from state X i ) gives the transition rates connecting all states X i with all states X j (q i j = 0 when i = j). To preserve the Markovian structure of the process, all the coefficients q i j are considered to be time-invariant.
The continuous-time master equation associated with p and Q is [32][33][34]: The Markov chain representing this process generates a sequence of state transitions X i → X j with probabilities Time reversibility of Eq. (1) results in the detailed balance condition: q i j p j = q ji p i (in matrix form, Qp T = Q T p), which allow us to express the above equation in matrix form as: where R i j = r i δ i j is a diagonal matrix (as a matter of notation, here summation is explicitly indicated by the symbol ' ' and the corresponding indices). The residence time in state X i , δt i can be shown to follow an exponential distribution with mean 1/r i [26,35], which allows us to define a probability Pr (δt i > t|ξ 1 (0) = X i ) = r i exp{−r i δt i } from which to sample the residence time using a uniform random variable ξ 2 : δt i = r −1 i log(1 − ξ 2 ). Solving the master equation (1) (or 2) fully characterizes the evolution of the system in time. However, in most applications of interest the total number N of states accessible to the system (dimension of the Markov space) is not known, resulting in undetermined Q matrices and thus precluding the use of analytical or deterministic methods to solve the master equation. In such cases, the problem is characterized by a discrete Markov chain with infinite dimension, and approximate methods must be used. When the following two conditions are met, discrete event methods can be used to solve the master equation: (i) the probability of each event depends only on the state reached in the previous event (i.e., the system is defined by a 1st order Markov process), and (ii) the system remains unchanged during a transition between two states (i.e., the time spent during a transition is negligible compared to the residence time). This is where the utility of the residence-time algorithm (also known as the kinetic Monte Carlo (kMC) method) resides. In particular, in rejection-free kMC (also known as 'BKL' after its original proponents [30]), once a transition between two states n → m is selected with the corresponding probability q nm / m q nm , it is executed and time is advanced by an amount equal to either 1/ m q nm or obtained by sampling a uniform random number χ ∈ (0, 1] as δt n = − log(1 − χ)/ m q nm . In this fashion, the system evolves in discrete time increments as the state space {X } is sampled, leading to a significant computational gain over standard Monte Carlo methods. The reader is referred to recent reviews for more information on the topic [34,36].

Kinematics
For a deformable body occupying a volume Ω 0 bounded by a surface ∂Ω 0 , a one-to-one mapping x (X, t) is assumed to exist between the position of material points in their reference position X and their current position x. The deformation gradient of this mapping, F = ∂x/∂X is typically decomposed multiplicatively into plastic and elastic contributions, F P and F E as [37]: where u is the displacement vector. The rate of change of F can be written as: where L is the velocity gradient. L additively decomposes into: where L E and L P are the elastic and plastic velocity gradients, respectively. In the small deformation limit (linearized kinematics), where H E and H P are the elastic and plastic distortions, i.e., L E =Ḣ E and L P =Ḣ P . If the only mechanism of plastic deformation at the crystal level is dislocation slip along specific crystallographic directions, then we can write [38]: whereγ α denotes the plastic slip rate of slip system α, and s α and n α represent the corresponding slip and plane normal directions (expressed in the original frame of reference). The tensor (s α ⊗ n α ) is known as the Schmid tensor. Under linear elasticity, the Cauchy stress can be obtained as a function of the elastic strain as, where C is the stiffness tensor of the crystal. For isotropic elastic systems, this second-order tensor can be written as a symmetric matrix C := (κ − 2μ /3) I ⊗ I +2μ1, where κ and μ are the bulk and shear modulus, and I and 1 are the second and fourth order identity matrices. The symmetricity of σ and C renders H E symmetric as well, which going forward we represent as ε E . Assuming an instantaneous elastic response, the rate form of Eq. (8) is: The system evolves in time subjected to the kinematic boundary condition: which combined with Eq. (9) results in the following relation, whereε 0 is a prescribed strain rate tensor.

Flow rule
The present CP model is rate dependent through the definition ofε P in Eq. (7). The shear ratesγ α are computed according to Orowan's equation: where ρ α , b α and v α are the dislocation density, Burger's vector modulus, and velocity in slip system α. Here, we consider bcc metals with only 1 /2 111 screw dislocations, i.e., b is unique and equal to a 0 √ 3/2, with a 0 the lattice parameter. In this work we are interested in low and intermediate homologous temperatures (< 0.3T m , where T m is the melting point), where thermally activated screw dislocation motion is represented by the following expression [39,40]: where ν 0 , h, λ α , w, and ΔH 0 are, respectively, the attempt frequency, the distance between Peierls valleys in the bcc lattice, the mean dislocation segment length for slip system α, the kink-pair width, and the kink-pair formation energy. Δτ α is the excess stress, obtained as the difference between the resolved shear stress (RSS) and the slip resistance, g α . The RSS is obtained from the Cauchy stress as: which is also known as the Schmid stress 3 , while g α includes all athermal sources of stress opposing dislocation motion. σ P is the Peierls stress. The time evolution of the dislocation density is modeled using a rate equation of the type: where ρ α mult and ρ α ann represent the net sources and sinks of dislocation line length. Each of these terms can be defined using a standard Kocks-Mecking model [39,42]: where it has been assumed that dislocation annihilation occurs spontaneously when dipoles meet within a critical spacing d edge ≈ b. λ α is the available dislocation segment length, which is generally expressed as: where d g is the grain size and ρ α f is the forest dislocation density, which is obtained for each slip system α from contributions of all other systems β [43]: The constitutive nature of ρ α comes via its dependence oṅ γ α , which depends on the stress through the dislocation velocity. This is an implicit dependence that will have numerical implications as will be discussed below. To close the model, an expression for g α must be provided. In the present case, we consider self and latent hardening only, both embodied in a term τ α h , i.e., g α ≡ τ α h , where: The first term inside the square root on the r.h.s. of the above equation represents the forest hardening, while the second one reflects the amount of self-hardening. ξ αβ are the elements of the hardening coefficient matrix, ξ [39,44,45], whose values, as well as the glide systems considered in the work, are given in Appendices A and B .

Casting the crystal plasticity model as a stochastic process
The Markov chain representing the process described by the master equation (1), (2) evolves the system forward in time by following a sequence of transitions whose rates depend solely on the initial and final states. As such, Eq. (1) is akin to an explicit time discretization of the crystal plasticity model, in which the properties of the system at time t + δt are obtained solely as a function of those at time t. With this, Eq. (9) can be recast as: The constitutive framework for the explicit form of the model remains unchanged. Additionally, the unknowns of the crystal plasticity problem remain the same. However, instead of solving iteratively for the ensemble of slip rates, the shear rates are determined in a sequential manner. The evaluation at time t +δt is based on the calculation of the available glide stress τ α RSS and the dislocation velocities at time t [46][47][48]. By analogy with (2), the above equation can be trivially written as a stochastic equation: In essence, this equation yields the time rate of the probability of finding the system under a given stress state. In the r.h.s. of the equation,ε 0 can be regarded as an external source term that is independent of the current stress state, whileḢ P acts as a stress absorber through conversion to plastic deformation. From Eqs. (11) and (20), Eq. (21) can be expanded to: While this equation is not strictly a master equation, i.e., one of the typeṗ = Q p, as in Eq. (2), a set of event rates can be extracted from Eq. (22) that govern the evolution of the elastic strain in time. The total rate to 'exit' a given stress Fig. 1 Sampling array of event rates in the SCP algorithm. The total rate r is partitioned between the prescribed applied strain rateε 0 and the total plastic distortion H P . In turn, H P is composed of N independent shear rates γ α . As such, there are always N + 1 different possible event rates to sample from, contained in an array whose index k runs from 0 to N . Note that, in general, the value of each shear rate γ α may be different to the rest of the shear rates during each time step. The dashed section refers to the possibility of using null events, discussed in Sect. 4.1 state defined by a total elastic strain ε E at time t n can be written as: and N is the number of independent slip systems.
In rejection-free kMC, the next event to be executed is selected as the k th process that satisfies: where k = 0, . . . , N . If the event selected corresponds to the rate q 0 , further sampling is carried out among the different nonzero components ofε 0 to determine which deformation to impose. The total time is then advanced as: where ξ 1 and ξ 2 are uniform random numbers in (0, 1]. Once an event is executed, the elastic strain is updated and, with it, the stress state. From the new stress state, resolved shear stresses are calculated, from which in turn plastic slip rates and dislocation densities are updated, thus closing the cycle. A schematic diagram showing the sampling of event rates and the total rate is given in Fig. 1. A complete algorithm suitable for solving Eq. (22) based on the above procedure is provided next.

A residence-time algorithm for elasto-viscoplastic CP problems
The crystal is defined by its dislocation density content, ρ α , in each of its N slip systems α and the grain size d g . Deformation conditions are given by the applied strain rate (tensor) ε 0 , the temperature T and the loading orientation o (generally aligned with a given crystal axis corresponding to the z direction). Thus our simulations are carried out with only one nonzero component of the imposed strain rate tensor, such that q 0 ≡ (ε 0 ) zz . The algorithm proceeds as a standard crystal plasticity algorithm until the calculation of the slip rates is completed. The stochastic part of the algorithm that serves as integrator is captured in lines 30 to 32. For comparison, the equivalent deterministic (explicit) algorithm is given in algorithm 2 in "Appendix C".

Physical bounds on problem time scale
Calculate: stress (tensor) increment Δσ = C : Δε E 4: Update: stress tensor σ = σ + Δσ 5: Initialize: q 0 =ε 0 6: Get: random numbers ξ 1 , ξ 2 ∈ (0, 1] 7: for α = 1, N do 8: Calculate: modulus: b = b α 9: Get: slip direction s α = b −1 b α , plane normal n α 10: Calculate: resolved shear stress τ α RSS = s α · σ · n α 11: Initialize: r t =ε 0 12: for β = 1, N do 13: Calculate: forest dislocation density ρ α f = ρ α f +ρ β |s β ·n α |. 14: Calculate: forest dislocation hardening g α = g α + ξ αβ ρ β |s β · n α | 15: end for 16: which are thus reached in much shorter times than those indicated in Fig. 2 . One way of establishing meaningful boundaries for δt may be by considering the physical validity range of the present algorithm, which is set by the condition (12)), i.e.: Taking time increments, the above expression becomes (σ P is constant): Implicit in this relation there is a time step condition that must be satisfied such that the time scale of the problem is reflective of the physical processes that participate in the evolution of the deformation of the material. As shown in Alg. 1, τ α RSS displays a linear dependence on δt, while τ h scales as O(δt 1/2 ) (through its dependence on the term β ρ β , β = 1 . . . N ). As such, one can rewrite the above inequality generically as: which results in the condition: where C 1 and C 2 are constants representing, respectively, the (elastic) stress rate, C 1 ≈ E ε 0 −γ p and the forest hardening rate C 2 ≈ μ γ p , whereγ p is again a generic plastic shear rate, i.e.: where E = 9κμ 3κ+μ is the Young's modulus. For simplicity, it is assumed that the shear rates are well captured by the term ρ 0 bv, where v is a temperature dependent dislocation velocity, set by Eq. (12). With this, the above expression becomes: This expression is fundamentally equivalent to that derived by Van der Giessen et al. [49] using a slightly different constitutive model. In the present model, for the 300-to-1000-K temperature range, v varies between 10 −7 and 10 −3 m·s −1 . By way of example, using material constants from Table 1, this results in critical time steps δt * between 0.05 to 100 s whenε 0 = 10 −3 s −1 .
As such, the physical time scale symbolized by δt * and the results from Fig. 2 (defined by the value 1/r n ) must be reconciled. While there are rigorous ways to connect both The impact of the choice of Δε * will be evaluated in the next section.

Verification of algorithm capabilities
In this section, we verify that the stochastic solver in Algorithm 1 is capable of reproducing the results of deterministic CP calculations under generic loading and temperature conditions. We first consider uniaxial loading along the [001] crystal orientation at 500 K and 10 −3 s −1 .
[001] loading results in eight active slip systems each with a Schmid factor of 0.409 and four inactive ones. As such, it is representative of multi-slip conditions conducive to latent hardening. In this work we consider single crystal deformation only, i.e., 1/d g ≈ 0. Values for the material parameters used in the calculations are given in Table 1. Figure 3 shows stress-strain curves for several Δε * . The shaded area represents the reference deterministic result (forward Euler method with dt = 0.1 s) for comparison. A prescribed strain rate of (ε 0 ) zz ≡ε 0 = 10 −3 s −1 was used. As the figure shows, the value of Δε * determines the magnitude of the stochastic oscillations, but not its steady state value, which is seen to consistently reproduce the deterministic solution. In particular, the hardening modulus -assumed to remain constant during plastic flow-is seen to be virtually identical in all cases, 0.96 GPa, irrespective of the value of Δε * . The average yield strength,σ Y , obtained from all the different Δε * cases is seen to slightly exceed the deterministic value, 1.9 vs. 1.8 GPa, respectively. The standard deviation ofσ Y scales linearly with Δε * , with a proportionality constant on the order of the Young's modulus E. This is consistent with having elastic strain increments that may 'overshoot' the true (deterministic) yield stress, which is then immediately followed by a plastic event. Figure 4 shows the evolution of the total dislocation density with strain corresponding to the cases shown in Fig. 3. The values shown in the curves are normalized to the initial dislocation density value ρ 0 (cf. Table 1). The dislocation density is seen to evolve in bursts of a magnitude commensurate with the value of Δε * . Although a small drift at high strains can be appreciated for the two cases with the largest scaling factor, the steady state value of ρ tot is in excellent agreement with the solution furnished by the deterministic solver. ρ tot is an integrated measure of the overall performance of the algorithm during plastic deformation. However, it is also of interest to look at the evolution at the level of the slip systems. Figure 5 shows results for the average shear rates under [001] axial loading at 500 K, defined as:  shown represent the standard deviation associated with the stochastic sampling. As in Figs. 3 and 4 , the deterministic solution is shown in the background as a shaded area. Clearly, the values associated with larger Δε * show much larger average plastic shear rates and more pronounced deviations from the mean value.

Exploring the parametric space
For generalization purposes, next we apply the method changing the three external variables in the simulations, i.e., temperature, strain rate, and loading orientation. Figure 6 shows a set of stress-strain curves at temperatures ranging from 300 to 800 K atε 0 = 10 −2 s −1 under [001] loading. All curves are for Δε * = 5.0 × 10 −5 . The temperature dependence of the yield strength (obtained as the 0.2%-offset stress) is given in the inset to the figure. Figure 7 shows the stress-strain response of the system as a function ofε 0 at 500 K under [001] loading. All curves are for Δε * = 5.0 × 10 −5 . The inset shows the strain-rate sensitivity (SRS) of the yield strength (obtained as in Fig. 6

Natural evolution of slip under heterogeneous conditions
The above subsections unequivocally show that the stochastic crystal plasticity (SCP) algorithm is capable of solving the same mathematical problem as standard crystal plasticity. However, as indicated in the introduction, this is not the main point of SCP. Stochasticity represents an intrinsic numeri- To showcase this feature of the model, here we study again uniaxial tests at 500 K, 10 −2 s −1 , and several orientations using the SCP method sampling one single slip event per time step. Our curves are compared to results obtained using the CP explicit solver by Kuchnicki et al. [46], which is also designed to induce asymmetric plastic flow by using a sequential cycle that prioritizes slip in systems with the highest excess stress (i.e., what the quantity Δτ α used in Eq. (12) represents here).
To demonstrate slip anisotropy, we take the [101] direction as representative of multislip conditions (four slip systems with a Schmid factor of 0.408 and eight inactive slip systems) to compare the deterministic approach by Kuchnicki et al. and the SCP method. In Fig. 10 we plot the plastic shear rates in all 12 slip systems for the [101] loading case shown in Fig. 9. As shown, the plastic slip that leads to the stress-strain curves for the [101] loading orientation in Fig. 9 is carried by a single slip system (α = 7 in Table 2) in the SCP case, and by two of them (α = 6, 9) in the deterministic case. We have confirmed that in each SCP simulation the probability that any one of the four active slip systems is activated is indeed the same (25%).

Computational performance
In general terms, the efficiency of the residence time algorithm is tied to the event search encoded in Eq. (24). When the dimension, n, of the sampling array is large, binary searches with a O(log(n)) nominal overhead are more efficient than simple O(n) linear searches. However, this is not a factor in the present calculations, where, at most, n = N + 2 with N = 12. Instead, the computational overhead may be defined as the CPU time invested in advancing a CP simulation by some amount of strain. Figure 11a shows the CPU cost per 1% strain 4 for the implementation of Alg. 1 on a 1.4GHz Intel Core i5 processor tested as a function of Δε * for a number of SCP simulations at 500 K under [101] loading at 10 −3 s −1 . As the figure clearly shows, an inverse correlation between the CPU time and Δε * is found. For the results shown in the figure, t CPU = 3.43 × 10 −5 x −0.97 [s per 1% strain]. This is not unexpected, since, as we showed above, there is a direct equivalence between Δε * and the effective time step used in the simulations. Using such correlations can be helpful in estimating the a priori CPU overhead of a SCP simulation.
For its part, Fig. 11b gives the same metric as a function of applied strain rate for a fixed value of Δε * = 10 −4 . A sharp drop is observed betweenε 0 = 10 −3 and 10 −2 s −1 , although this is likely to be somewhat dependent on the value of Δε * chosen in each case.

Discussion
In essence, the present method provides a new approach for time integration of crystal plasticity models. The solver itself is the well-known residence time algorithm, which relies on stochastic sampling of a set of transition rates that determine the evolution of the system through a sequence of states with the correct probability. KMC methods require concurrent sampling of all event rates, and thus the method is suited for explicit time discretizations of the crystal plasticity equations (whenε 0 and the differentγ α are strictly independent of one another). However, this is more a technicality in the definition than a limitation in the computational sense, as the SCP approach is not constrained by the same stability considerations of standard (deterministic) methods.
The discrete nature of the method naturally leads to oscillations in the plastic response of the system. It is worth emphasizing that, while numerical in origin, these oscillations are physical and can be ultimately linked to the discreteness of dislocation slip at the lattice level, by way of integer Burgers vector amounts. The magnitude of the oscillations is encoded in a parameter Δε * , which here is user-defined under certain physical restrictions (Sect. 2.5). A way to quantitatively interpret these oscillations may be by considering Δε * to scale with the ratio b/L, where b and L are the Burgers vector's modulus and a generic specimen's dimension, respectively. In bulk materials, this ratio becomes Fig. 10 Normalized plastic shear rates as a function of plastic strain for symmetry-breaking slip evolution modes of the SCP algorithm and an explicit deterministic model [46] (a) (b) is seen in small-scale specimens such as free-standing nanopillars or cantilever nano-beams [50][51][52][53]. This, and other features of the model, including several potential advantages over standard explicit CP approaches, are further elaborated on below.

Physical time scale defined through 1" *
In Sect. 2.5, a parameter Δε * is introduced to restrict the time steps that would naturally emanate from Eq. (23), graphically represented in Fig. 2. However, a formal connection between Δε * and δt * , in Eq. (31), was not established at the time.
A rigorous way to demonstrate this connection can be by adding an extra rate to r n in Eq. (23) to extend the exit rate of the system with an event that does not alter its state at a given time. This 'null' event process, which has been applied successfully in the context of parallel kMC algorithms [28,54], is characterized by a rate r null that can be used at will to 'slow' down the evolution of the system while preserving the correct kinetics. Mathematically: which leads to the condition: r null > 1 δt * − r n . The pictorial representation of r null is shown in Fig. 1 as an extra entry (dashed gray line) in the array. It is trivial to show that the relationship among Δε * , r n , and r null is: Δε * = r n r n + r null or, alternatively, Δ * /r n < δt * , i.e., Δε * , r null , and δt * can be thought of as being interchangeable parameters. In any case, these relationships show that one can rigorously connect the mathematical representation of the system's kinetics to the physical time evolution under deformation.

Asymmetry in slip rates with identical orientation
It is well known that from a purely geometric point view, there is a redundancy in establishing the number of active slip systems that contribute to deformation under highly-symmetric loading conditions. Mathematically, this problem of slip indeterminacy is intrinsically found in rate-independent CP formulations, and can be partially addressed by, e.g., imposing additional constraints on the internal variables or the stress space of the system [13,55,56]. While the non-uniqueness of the solution is no longer an issue in viscoplastic (rate-dependent) crystal plasticity models, the plastic response still displays symmetries in slip that are not observed in real materials' deformation [57][58][59].
The SCP model presented here offers the potential to introduce variability in the solution in a natural way, thus intrinsically enabling the possibility for non-symmetric activation of slip systems according to just the Schmid criterion. This symmetry-breaking capability emanates directly from the solver, i.e., without the need to add it externally, which means that it can lead to inhomogeneous plastic deformation naturally. This inhomogeneity emerges from the inherent stochastic treatment of the plastic shear rates, and thus may be seen as having a direct correlation to real deformation situations. From a thermodynamic standpoint, several works have established a link between the heterogeneity of plastic deformation and the system's entropy production based on maximum-dissipation criteria [60,61]. Indeed, stochastic approaches have the potential to access an increased number of states. While beyond the scope of the present work, it would be of interest to augment the current model by introducing quantitative metrics of entropy production directly connected to the stochastic variability and the symmetrybreaking features of the SCP model.

Numerical stability and computational cost
While the computational overhead of the SCP method has been assessed in Sect. 2.2 and Fig. 11, any meaningful discussion of the CPU cost of the SCP method must necessarily involve a comparison with a standard explicit (forward Euler) integrator for CP problems in equivalent conditions. Figure 12 shows the computational cost (in seconds per 1% strain) of a reference crystal plasticity problem ([101] loading direction atε 0 = 10 −3 s −1 ) as a function of temperature simulated both with the SCP model using several values of Δε * and a standard deterministic CP model based on a forward Euler integrator. The numbers associated with each data point of the deterministic calculations (in blue) indicate the maximum dt attainable to maintain numerical stability 5 .
There are several remarkable features that emerge from the analysis presented in the figure. We start by noting that the SCP method is intrinsically stable under all conditions, a virtue of the residence-time algorithm, while the numerical stability of the deterministic solver is highly dependent on temperature and strain rate. Second, the CPU cost of the stochastic algorithm is practically independent of temperature. Third, and most important for the purposes of this discussion, the deterministic solver's numerical stability requires increasingly smaller time steps (labeled as 'δt max ' in the figure) as the temperature increases. For this reason, at 600 K we see a sharp upturn in the CPU cost, which drives its computational overhead far above that of the SCP model. The fact that the algorithm samples one event per time step helps stabilize the numerical solution and allows for longer time steps at high temperature compared to a standard forward Euler approach. While it is common practice to add modifications to the Euler method to improve its stability [9,49], we consider this a very attractive potential advantage of SCP over standard explicit solvers.
We close Sect. 4.2 reiterating the potential advantages of SCP over deterministic (both explicit and implicit) approaches. One obvious advantage of using our method at all temperatures (including room temperature and below) is the added benefit of capturing plastic fluctuations naturally. Thus, in cases where those are of interest (cf. Sect. 4.2.1), our approach is useful from a physical point of view. Numerically, at low temperatures our method is consistent with the behavior of kinetic Monte Carlo models versus deterministic solvers in other fields of physics [62][63][64]. Rather than being attributable to stochastic solvers being slow, the advantages of deterministic solvers have to do more with their ability to be stable using large time steps. In general, the consensus is that for simple models (with low physical complexity) deterministic models will always outperform stochastic ones. However, as the physical model (so-called 'material model' in crystal plasticity) grows in complexity, stochastic solvers become relatively more efficient. Note that such complexity may be introduced also via the external conditions, as is the case in this work, when thermally activated mechanisms become dominant at higher temperatures.

Conclusions
We conclude this paper with a list of our most important findings: -We have developed a stochastic solver for crystal plasticity models based on the residence-time algorithm. The method strictly works for explicit problems, when the total strain rate and the individual plastic shear rates can be considered independent from one another.
-The SCP model is intrinsically numerically stable without the need of any extra procedures. Changes in the values of the external variables manifest themselves in terms of the magnitude of the oscillations of the solution, not on its stability. -The SCP model naturally breaks the symmetry of plastic slip by sampling among the active plastic shear rates with the correct probability. This can lead to phenomena such as plastic localization without needing to add any ad hoc treatments to the model. -All variables kept the same, the computational overhead of the SCP method scales inversely with Δε * (as ∼ 1/Δε * ) andε 0 . The CPU cost is insensitive to temperature for the crystal plasticity model employed here. -For a fixed prescribed total strain rate and loading direction, the SCP model becomes more efficient than a standard forward Euler approach at T > 600 K. It is expected that similar transitions exist for other crystal orientations and strain rates.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A Crystallographic slip systems for bcc lattices B Latent hardening coefficients
See Tables 3 and 4.

C Explicit deterministic algorithm
The reference deterministic solver for the explicit crystal plasticity model used here is given below. It is based on a simple forward Euler method without any modifications to improve its convergence behavior (first order accuracy, O(δt 2 )).