Stochastic Crystal Plasticity Models with Internal Variables: Application to Slip Channel Formation in Irradiated Metals

Stochastic crystal plasticity models are originally introduced to study slip avalanche phenomena that are ubiquitous features of the microscale plasticity of crystalline solids. Herein, a method is proposed to couple such models to the evolution of internal variables to account for microstructural hardening and softening phenomena. Specifically, strain hardening is described in terms of a Kocks–Mecking‐type dislocation density and the structural softening of irradiated metals in terms of the density of irradiation‐induced point defect agglomerates, which are cut and eliminated by moving dislocations. The interplay of both effects results in the formation of defect‐free slip channels. Critical conditions for slip channel formation are formulated, the statistical morphology of the ensuing slip channel patterns is investigated and compared with surface observations. Finally, the magnitude and nature of stress concentrations that emerge if slip channels interact with platelet‐like hard inclusions are discussed.


Introduction
Crystal plasticity (CP) governed by the collective motion of dislocations is, on the microstructural scale, characterized by very significant fluctuation phenomena as manifested in macroscopic samples by the emergence of intermittent acoustic emission signals, attributed to dislocation avalanches with scale-free statistics, [1,2] and in micrometer-scale samples by directly observable strain bursts. [3,4] To computationally model these phenomena, two alternative approaches have been pursued. On the one hand, discrete dislocation dynamics (DDD) simulations have been performed in two [5] and three [4] dimensions and the dislocation avalanches observed in these simulations have been statistically analyzed. On the other hand, stochastic CP models have been introduced, [6,7] which generalize standard CP to account for microstructural heterogeneity and disorder and for the ensuing intermittency and jerkiness of plastic flow. These models have until now been only used to describe intrinsic fluctuations of the plastic shear strain rates and the resulting spatiotemporal fluctuations of strain; in this sense, they are akin to standard CP models which consider plastic shear strains on active slip system(s) and their rates, but no other internal variables. However, many CP models couple a description of plastic shear strain evolution with the evolution of internal variables describing the material microstructure, notably dislocation densities (see, e.g., the overview of Roters et al. [7] ). In this manner, it is possible to account for plasticityinduced microstructure evolution and formulate physically based models for strain hardening. In the present study, we couple for the first time a stochastic CP model with differential equations describing the evolution of defect and dislocation densities, which thus assume the character of stochastic differential equations. This allows us to describe the formation of slip channels in irradiated metals and model and analyze the ensuing, strongly heterogeneous slip patterns.
Slip channels [8] emerge during plastic deformation of irradiated crystals containing high densities of irradiation-induced defect agglomerates such as voids, bubbles, [9] or stacking-fault tetrahedral. [10] The channels form thin lamellar regions parallel to the primary slip plane and are associated with a very significant localization of slip on the microstructure level. [8] As a consequence, in locations where slip channels terminate at grain boundaries such that the localized slip is not (or not fully) accommodated by slip in the adjacent grain, large stress concentrations may emerge and lead to grain boundary fracture. [11] This effect is believed to be a key factor in irradiation-assisted stress corrosion cracking and thus in irradiation embrittlement. [12] On the microscopic level, slip channels arise because moving dislocations can absorb point defect agglomerates of both vacancy and interstitial type. The corresponding dislocation processes have been investigated in a number of studies, using both in situ experimentation (see, e.g., previous studies [13,14] ) and molecular dynamics (MD) simulation (see, e.g., previous studies [14,15] ). Since irradiation-induced defect agglomerates act as dislocation obstacles, their removal leads to local softening, which facilitates dislocation motion and slip localization in embryonic channels. In fully developed slip channels, this process is counteracted by dislocation multiplication and dislocation-dislocation interactions, leading to mutual trapping of dislocations, e.g., in the form of dipolar and multipolar bundles and mats. [8] The ensuing strain hardening delimits the channel slip and is supposed to control the channel width. [8] From a modeling point of view, the elementary processes of dislocation-irradiation defect interaction have been studied extensively using MD simulations. [14,15] In multiscale frameworks, local rules for obstacle-dislocation interactions have been parameterized on the basis of MD results or in situ observations and implemented into DDD simulations. [16,17] An alternative approach, which has the advantage of giving access to much larger system sizes (and thus polycrystal phenomena), consists of CP simulations which use dislocation densities and irradiation defect densities as internal variables. [18][19][20][21] Importantly, to correctly reproduce the strongly localized slip distribution associated with channel slip, such approaches have to introduce some kind of microstructural randomness into either initial conditions or evolution equations. In a previous study, [18] this was done in terms of random fluctuations of the initial density of irradiation obstacles. Here, we follow the same general line of modeling but consider randomness in terms of the well-documented intrinsic fluctuations that characterize dislocation flow even in the absence of irradiation hardening. The article is organized as follows. In Section 2, we first briefly discuss the general approach of stochastic CP simulation and its coupling with evolution equations for density evolution of dislocations and irradiation defects. Issues of numerical implementations are also discussed. Section 3 presents simulation results, first for free channel formation and second for slip channels that are blocked by extended obstacles. We evaluate the ensuing stress concentrations and show that the blocking of slip channels leads to stress and slip patterns that may be envisaged as shear anticracks. We conclude with a discussion and outlook in Section 4.

CP with Internal Variables
CP envisages plastic deformation in terms of shear strains on crystallographic slip systems. For simplicity and conceptual clarity, we use a small strain formulation. A slip system ϕ is then defined by a fixed crystallographic unit vector s ϕ defining the direction of shear parallel to a lattice plane with unit normal n ϕ . The plastic strain tensor ε pl is expressed in terms of projection tensors and (scalar) shear strains γ ϕ on the active slip systems as The thermodynamic driving forces work to conjugate to the shear strains γ ϕ and give rise to the resolved shear stresses The stress tensor σ fulfills the constitutive equation with C being the tensor of elastic constants and is determined for a given plastic strain field by solving a standard elastic-plastic boundary value problem. The CP framework is complemented by constitutive equations, which relate the shear strain ratesγ ϕ on the active slip systems to the corresponding driving forces τ ϕ , e.g., where τ k c ðfγ j gÞ are flow stresses on the different slip systems. The rate-independent limit is recovered by considering large values of the rate exponent m. Strain hardening can be expressed either by directly specifying the dependency of the flow stresses on the shear strains, e.g., in form of a hardening matrix, or by relating the flow stresses to internal variables that in turn evolve with strain.
In internal-variable approaches to CP, one considers physically motivated hardening laws which relate flow stresses to the densities of defects, such as dislocation densities on the different slip systems or densities of irradiation-induced obstacles, which in turn evolve with strain. As an example, we consider a Taylor-type hardening law which relates the slip system-specific flow stresses to slip system-specific dislocation densities ρ ϕ where α % 0.25 is a nondimensional factor (Taylor factor), μ the shear modulus, b the length of the Burgers vector, and h ϕψ is a latent hardening matrix. In a dislocation density-based formulation, the rate law, Equation (5), may be modified using Orowan's relatioṅ To yield a complete description of plastic flow, the constitutive Equation (6) and (7) must be complemented by equations for the strain evolution of the internal variables (dislocation densities). These may, for instance, be taken of the Kocks-Mecking type In this equation, β À1 is the dislocation mean free path in units of dislocation spacings, and y ρ is an effective cross-section for dislocation annihilation. If two dislocations of opposite sign (n 3D: of opposite line orientation) meet on parallel slip planes at a slip plane distance below y ρ , then they are supposed to annihilate by cross-slip or climb. For edge dislocations in irradiated Cu, Essmann and Mughrabi [22] give a value of y ρ ¼ 1.6 nm, though depending on cross-slip propensity, the value for dislocations with near-screw orientation may be much higher.

Stochastic CP
Stochastic CP retains the general CP framework but differs in the interpretation and evaluation of the flow stress and shear strain rate. Whereas standard CP considers these quantities as continuous differentiable fields, which obey deterministic equations of evolution, stochastic CP adopts an essentially discrete and stochastic formulation. The underlying physical ideas are as follows: 1) Dislocation plasticity is controlled by the formation of localized, metastable configurations, with spatial extensions of the order of the mean dislocation spacing. In single slip, these are dipoles and multipoles; in multiple slip, these are dislocation junctions. 2) At zero temperature, dislocations are released from such configurations once the local shear stress exceeds a critical threshold. We are thus dealing with threshold dynamics. 3) After a threshold is overcome, dislocations rapidly sweep a finite area with linear dimensions of the order of the dislocation spacing, producing a finite strain increment. 4) As the microstructural information available on the CP scale is necessarily incomplete, statistical uncertainty about the actual microprocesses requires a statistical description of threshold stresses and strain increments.
To implement these ideas into a CP setting, stochastic CP models as proposed in previous studies [6,7] adopt a discrete, automaton-type dynamics. Space is discretized into grid units ("mesoscopic elements") with an extension of one dislocation spacing. Each grid unit i carries a local plastic strain, which is evaluated from shear strains on the different slip systems as in standard CP; hence, Equation (1)-(4) of the CP framework is retained. However, modifications are applied in the following respects:

Statistical Distribution of Flow Stresses
The critical stressτ ϕ i required to move a dislocation of slip system ϕ across a mesoscopic volume element i is controlled by the configuration of dislocations within this volume, in particular by the spacing and arrangement of forest dislocations, with which a moving dislocation necessarily forms junctions. The stress required to unzip a dislocation junction is inversely proportional to the junction length, hence to the characteristic dislocation spacing. This is the origin of the Taylor relationship in multiple slip. (NB: We note that the argument can be transferred to single slip, if one replaces "junctions" by "dipoles" and "junction length" by "dipole height." For generic mathematical arguments underpinning the scaling relations used here, see previous study. [23] ) However, junction lengths are not uniform but, as demonstrated by DDD simulations, [24] they have a very wide statistical distribution. Now the key problem is that in a CP simulation which operates on scales above the resolution of individual dislocations, we cannot know the actual junction length in a given grid element, and hence its strength. We only know the average strengths τ ϕ c of elements (the average is here understood as an average over many statistically equivalent elements with the same dislocation densities), which obey the Taylor relationship Equation (6). The actual local strength values have to be inferred statistically, accounting for the fact that they will deviate from the average values. We therefore have to develop a rule for such statistical inference that allows us to assign local flow thresholds to our grid elements.
What we know is the mean junction length, which is proportional to the dislocation spacing, and thus the mean strength, which must obey the Taylor relationship. Based upon the mean strength τ ϕ c , we can estimate the probability distribution Pðτ ϕ i jτ ϕ c Þ of local strengths using an information theoretical argument [25] (see A1, Supporting Information) as Then, instead of assuming that each element has the same threshold, we assign to them independent random strength values from the distribution Pðτ ϕ i jτ ϕ c Þ(see A2, Supporting Information).

Evolution of Strain
On the microscale and in the rate-independent (low-strain-rate) regime, dislocation motion does not proceed smoothly. Rather, dislocations are trapped in metastable configurations, in which they remain until the local stress reaches a threshold level. They then move rapidly until they are trapped into a new metastable configuration, such that the "flight" time between obstacles is negligible as compared with the time spent in obstacle configurations. We represent this dynamics by estimating the strain produced locally when a dislocation moves between two trapping configurations and assuming that this strain is produced instantly once a local threshold is overcome. We thus replace Equation (7) by the discrete incremental rule where Φ i is a random number of unit mean, with local stressdependent corrections to ensure positive definiteness of the plastic work. For motivation of this rule, see A3, Supporting Information.

Evolution of Internal Variables and Local Flow Thresholds
The local flow thresholds are functions of local strain. This concerns both their mean values τ ϕ c and the actual threshold in a given element i. The mean thresholds are evaluated according to Equation (6), or analogous equations, from internal variables whose evolution is still given by Equation (8). However, since strain changes according to Equation (10) are finite, the same is true for dislocation density increments. Once local dislocation density increments are evaluated subsequent to a local strain increment, we recompute the mean threshold and then assign a new stochastic threshold according to the distribution in Equation (9).

Simulation Protocol
A typical simulation protocol then reads as follows. 1) Once the computational grid is defined, local threshold values are assigned to all slip systems on all sites in the form of independent random variables distributed according to Equation (9). 2) The external load (imposed either in form of boundary tractions or displacements) is increased until, at one site, the local resolved shear www.advancedsciencenews.com www.aem-journal.com stress on one or more slip systems (the active slip systems) exceeds the local, slip system-specific threshold stress. Such a site is called active.
3) The local shear strains on the active slip systems on an active site are increased by finite amounts according to Equation (10), the local internal variables and mean flow stresses are recomputed according to Equation (6) and (8), and a new threshold value is assigned to the active site from the distribution in Equation (9). 4) The stress field in the entire system is recomputed. 5) As a consequence of the local stress changes, further grid sites may become active. Strain increments are applied on the active slip systems on these active sites, local thresholds on the active sites as well as global stresses are recomputed, and the procedure is repeated until all sites are inactive. 6) The external load is again increased by the minimum amount required to activate at least one site. 7) The simulation is terminated once a prescribed end stress or end strain has been reached. The ensuing dynamics proceeds in avalanches during which multiple sites become active. Consistent with experiment, these avalanches exhibit power law size distributions [2,6,7] which are, at large avalanche sizes, delimited by a cutoff that depends on system size and on the strain hardening coefficient, [26] in agreement with experimental observations and DDD simulations. [4] The resulting slip distributions imply a self-affine surface morphology. [2,6,7] Thus, stochastic CP captures essentials of the fluctuation phenomena associated with dislocation slip. We now adopt this framework to describe the formation of slip channels in irradiated metals.

Model and Computational Implementation
In our model of slip channel formation, we consider the earlystage deformation of a fcc single crystal where only a single slip system is active. Internal variables considered are dislocation density ρ (we drop the slip system-specific label) and volume density of irradiation-induced defect agglomerates η. Agglomerates act as strong obstacles that pin dislocations individually, provided the slip plane is within an effective distance d from the obstacle centre. The mean flow stress is then, in a grid element with the internal variables (ρ,η), evaluated as The equation of evolution for the dislocation density is of Kocks-Mecking type, and the evolution of the obstacle density is described by a simple first-order decay of the obstacle concentration with local strain In Equation (13), it has been assumed that an obstacle is annihilated, swept, or absorbed by a moving dislocation if it is located within a distance of less that y η /2 from the dislocation slip plane. From the local flow stresses, statistically distributed local flow thresholds are evaluated as detailed in Section 2 and A2, Supporting Information. We implement the model on a regular, 2D simple cubic computational grid of size L Â L. Periodic boundary conditions (BC) are used in both x and in y directions, and loading is conducted by imposing an affine shear ε xy ¼ ε yx ¼ γ ext on the periodic supercell. The slip direction of the only active slip system is identified with the y direction and the slip plane normal is in x direction; hence, we consider a plane strain situation. The system develops, due to the statistically distributed flow thresholds, heterogeneous plastic strain patterns. The corresponding internal stresses are calculated using a spectral approach, i.e., we convolute the plastic shear strain field with a Green's function as detailed in previous studies. [6,7] In our simulations, we assume material parameters (shear modulus μ, Poisson number ν, and Burgers vector length b) of copper. The Taylor parameter α is taken to be 0.25, and the initial dislocation density is taken to be ρ 0 ¼ 10 12 m À2 , which gives for a crystal without irradiation-induced obstacles a low initial flow stress (critical resolved shear stress) of 2.5 MPa. Our default grid size is L ¼ 512. With the understanding that one grid cell has a linear dimension of ρ 0 À1/2 ¼ 1 μm, a physical system size of about 0.5 mm 2 is implied. The dislocation multiplication parameter is chosen as β ¼ 0.02 to yield (in the absence of irradiation-induced obstacles) a hardening coefficient of μ/400 as typical of the end of hardening stage I in fcc metals. The dislocation annihilation distance is taken to be y ρ ¼ 1.6 nm, as determined by Essmann and Mughrabi. [22] Regarding irradiationinduced obstacles, we consider stacking-fault tetrahedra with a size of 50 nm, which we identify with the parameter d, and set the annihilation cross-section to y η ¼ d/2 ¼ 25 nm. We measure their density in units of η 0 ¼ α 2 ρ 0 /d which is the obstacle density for which according to Equation (11) the critical resolved shear stress is exactly twice as large as with the initial dislocation density alone. An overview of all parameters is given in Table 1.

Stress-Strain Curves and Microstructure Patterns
Stress-strain curves resulting from strain-controlled loading of a simulated volume are shown in Figure 1, together with corresponding obstacle density patterns recorded at 1% total shear strain. It is clearly shown that slip channel formation requires a minimum density of irradiation-induced obstacles to be present. In the present case, for η(0) ¼ η 0 , obstacles are removed homogeneously and the slip pattern does not differ appreciably from the unirradiated material. At η(0) ¼ 4η 0 , large fluctuations of the obstacle density are visible which, however, do not yet amount to fully developed slip channels. At higher densities (η(0) ¼ 9η 0 ), sharply localized channels emerge in which the initially present obstacles are almost completely removed. This picture then remains unchanged when the irradiation-induced obstacle density is further increased. The slip channel patterns are stochastic, without apparent regularity or periodicity. The question of correlations in the slip channel patterns will be discussed in the next subsection. Formation of slip channels is also visible on the stress-strain curves. Increasing obstacle density implies an increased flow stress and the emergence of a stress plateau at the beginning of the stress-strain curve (hardening stage I). The extension of the stress plateau is proportional to the initial irradiationinduced increase in the critical resolved shear stress. As shown in the inset of Figure 1, within the plateau, the local strength first decreases, then increases as a function of local strain. However, the system avoids the softening regime by jumping, as soon as yielding occurs locally, to the hardening branch of the stressstrain curve. Macroscopically, this sudden process of forming a slip channel is accompanied by a small load drop, and the stress-strain curves therefore assume a serrated shape in the regime of slip channel formation (see curves 3 and 4 in Figure 1). On the microstructure level, a kind of two-phase microstructure emerges where an essentially undeformed matrix with close to the initial density of irradiation obstacles coexists with slip channels where obstacles have been completely removed. The strain in a fully developed slip channel equals the plateau strain, and global deformation proceeds by increasing the slip channel fraction. Once the obstacles are removed everywhere, the plateau is followed by a hardening regime of constant slope (hardening stage II), which does not depend on the initial obstacle density.

Avalanches
Stochastic CP models were originally devised to model and analyze the avalanche dynamics, which is a characteristic feature of dislocation-mediated CP. It is therefore interesting to investigate how the microstructural instability associated with slip channel formation affects the avalanche dynamics. Avalanche size distributions for an obstacle density of η(0) ¼ 9η 0 are shown in Figure 2, where distributions of avalanches are plotted separately for hardening stages I and II.
We first focus on the hardening regime (hardening stage II). In this regime, the avalanche statistics does not depend on the initial irradiation defect density; it is in fact identical with the statistics for an unirradiated crystal. The statistics exhibits the same behavior as reported in previous works: the distribution follows a power law with an exponent of about 1.35 and an exponential cutoff. This cutoff increases in proportion with system size L-in fact, it can quantitatively be described by the relations Figure 1. Stress-strain curves, parameters as in Table 1, with respect to different initial defect densities from lowermost to uppermost curve: η ¼ 0, η ¼ η 0 , η ¼ 4η 0 , η ¼ 9η 0, η ¼ 16η 0, and η ¼ 25η 0 ; inset: stress-strain curve for η ¼ 9η 0 (red), together with local flow stress versus local strain (blue) as obtained from integration of Equation (11)- (13). Colorscale images show defect density distributions at 1% strain for different initial defect densities: 1) η ¼ η 0 , 2) η ¼ 4η 0 , 3) η ¼ 9η 0 , and 4) η ¼ 16η 0 , defect density on the colorscale is given in units of the respective initial density.
www.advancedsciencenews.com www.aem-journal.com determined by Zaiser and Nikitas [26] from a model without irradiation defects.
Interesting features arise, however, when the avalanche statistics in the plateau regime (hardening stage I) is analyzed. This is the regime when slip channels form, and their formation has direct consequences for the avalanche statistics: The avalanche size distribution is now no longer a scale-free power law but exhibits a pronounced peak at the size of the largest avalanches, which again increases in proportion with system size. In fact, this characteristic avalanche size directly relates to the characteristic plateau strain: A single slip channel has a width of one elementary cell and hence the 2D model comprises 1 Â L sites. For the parameters in Figure 2, the plateau strain of about 3% corresponds to about 120 elementary events (i.e., 120 dislocations passing each elementary cell); hence, formation of a slip line involves about L Â 120 elementary events, leading to characteristic avalanche sizes of S ¼ 7600, S ¼ 15 200, and S ¼ 30 400 for L ¼ 64, 128, and 256, respectively. A comparison with Figure 2, right, demonstrates that these numbers match well the location of the respective peaks in the avalanche size distributions. We may thus conclude that slip channel formation introduces a characteristic avalanche strain that changes the scale-free statistics normally associated with avalanche size distributions in CP.

Surface Morphology
Slip channel formation leads to a strongly heterogeneous surface strain distribution ("coarse slip") with huge surface steps and significant surface roughening. To assess the impact of slip channels on the surface in our model, which uses periodic BC, we assume that the strain distribution is the same in the bulk and on the surface. We can then obtain a proxy for the strain pattern in form of the displacement pattern on the side surface of our simulation cell, which we obtain by integrating the corresponding strain pattern. Thus, a 1D displacement profile is obtained as We analyze the stochastic morphology of the surface profile by determining the first-order structure function defined as where the integral is evaluated by making use of the periodic BC. An example of a surface strain pattern, the corresponding surface profile, and the corresponding height-height correlation function is shown in Figure 3. For self-affine surface profiles, the double-logarithmic plot of Φ(l) versus l shows a linear scaling regime, Φ ¼ Φ 0 (l/l 0 ) H , where l 0 is without loss of generality taken to be 1 μm, H is the Hurst exponent of the surface profile, and Φ 0 characterizes the absolute magnitude of surface height fluctuations. Figure 4 shows, for different initial obstacle concentrations, plots of H and Φ 0 versus strain (to obtain adequate statistics we average, for each data point, over typically ten profiles). We see two different behaviors: For unirradiated samples, the value Φ 0 is small and approximately strain independent, and the Hurst exponent fluctuates around values close to H ¼ 0.5. Such a value, typical of a random  Table 1, initial defect density η ¼ 9η 0 , different system sizes; left: avalanche size distributions in the hardening regime (strain >5%), right: avalanche size distributions in the slip channel regime (0.5% < strain < 3%). In irradiated samples where we find slip channels, in contrast, the characteristic fluctuation magnitude Φ 0 is much higher from the onset and then further increases until halfway into the stress plateau, up to a strain level where about half of the specimen volume is occupied by slip channels. The peak values of Φ 0 , which are up to two orders of magnitude above those in unirradiated material, increase with increasing density of irradiationinduced obstacles. From the peak value onward, the fluctuation magnitude drops and, at the end of the stress plateau, it becomes comparable with the values found in unirradiated crystals. Analogous behavior is observed for the Hurst exponent. In the slip channeling regime, initial Hurst exponents are high (H % 0.8) and they fluctuate around this level throughout the stress plateau regime, after which they drop to the significantly lower level H % 0.5 that is characteristic of unirradiated crystals.
Hurst exponents above 0.5 are indicative of long-range correlations in the slip pattern, [23] and we can thus conclude that slip channel formation not only leads to a significant surface roughening, which increases with irradiation, but also modifies the surface patterns and promotes long-range slip correlations.
Finally, a multiscaling analysis [27] of simulated surface profiles is included in A4, Supporting Information.

Interaction of Slip Channels with Extended Obstacles
The interaction between slip channels and extended dislocation obstacles such as grain boundaries is supposed to play an important role in irradiation embrittlement. Here, we study a simple, idealized obstacle configuration in the form of an extended, platelet-like inclusion which is oriented perpendicular to the slip plane. The inclusion is assumed to have the same elastic constants as the surrounding material but is assumed impenetrable to dislocations, hence plastically nondeformable. Figure 5 shows how such inclusions modify the slip channel patterns and internal stress patterns. In the regime of slip channel formation, the presence of a platelet like inclusion leads to a strong modification of the slip channel pattern: slip channels no longer form in a stochastic manner but emerge in close mutual proximity to form a shear band which gradually widens as the overall strain increases. Where the shear band intersects the inclusion, large internal stress concentrations emerge which cause a significant additional shear stress to act on the inclusion. On the edges of the intersection region, by contrast, the internal stress is negative, leading to a curious shear stress pattern, which resembles the stresses around a mode-II crack with a reversed sign.
The magnitude of the internal stress concentration on the inclusion can be evaluated by determining, as a function of the externally applied strain, the peak internal stress value. Results are plotted in Figure 6 for different densities of irradiationinduced obstacles. If the conditions for slip channel formation are not met, the peak internal stress value increases in approximate proportion with the external stress, following the hardening curve for a defect-free crystal. In the regime of slip channel formation, by contrast, the initial stage of deformation is   Figure 1 for the same parameters.
www.advancedsciencenews.com www.aem-journal.com characterized by significantly enhanced internal stress concentrations which increase with increasing density of irradiationinduced obstacles (and hence with increasing critical resolved shear stress) and are consistently much higher than the externally applied stress. The shear stress concentration is largest immediately after yield and then reaches a plateau before finally dropping, at the beginning of the linear hardening stage, to a level comparable of an unirradiated crystal. Note that such significant stress concentrations emerge despite the fact that, in the present simulations, plastic accommodation by slip on the opposite site of the obstacle is not impeded.

Conclusions
We have implemented a model of slip channel formation into the framework of stochastic CP models. We focused on the earlystage of deformation (hardening stage I) of fcc single crystals, considering irradiated Cu as a paradigmatic model material.
The qualitative predictions of the model match well the classical observations on irradiated metals: 1) Irradiation affects mainly the first stage of deformation, whereas subsequent strain hardening is little affected (see previous study [28] ); 2) Slip channel formation requires a critical density of irradiation obstacles, i.e., a critical irradiation dose (see previous study [29] ); and 3) Slip channels are associated with large surface steps and a coarse overall slip pattern. [28] Once slip channels form, it has significant consequences for external surface roughening and internal stress concentrations at obstacles. At the surface, the slip patterns show a transition from fine slip to coarse slip, resulting in an overall surface roughness that increases by more than a factor of 10 once slip channels form. The surface roughness becomes even larger when the density of irradiation obstacles further increases, such that global surface height fluctuations may exceed those in the absence of irradiation-induced obstacles by up to two orders of magnitude. Moreover, the slip patterns develop long-range correlations as manifested by a Hurst exponent H % 0.75 of the surface profiles that significantly exceed the value H ¼ 0.5 expected for the uncorrelated activation of slip. These predictions can be easily checked experimentally using standard surface characterization equipment. [23] If extended dislocation obstacles such as platelet-like inclusions are present, slip channels lead to significant internal stress concentrations at these obstacles which increase with increasing degree of irradiation but are absent as long as the critical conditions for slip channel formation are not met. Extended obstacles may lead to the correlated formation of slip channels, which further exacerbates the stress concentration.
Our mesoscale model allows the observation that slip channel formation is associated with coarse slip, surface roughness, and long-range correlations. Such phenomena live on length scales, which are not accessible with MD and DDD simulations. As slip channels have widths of the order of hundreds of nanometers, MD simulation can only describe the elementary processes of, e.g., dislocation motion, dislocation obstacle interaction, and obstacle removal annihilations. At the same time, DDD simulations are capable of modeling a limited number of slip channels, but not large-scale collective phenomena, involving multiple slip channels in interaction with each other as control, e.g., the surface morphology.
All the aforementioned phenomena are, in the present simulations, restricted to the first stage of deformation (hardening stage I) which is, in irradiated materials, characterized by little or no global hardening. Once this stage is over, the behavior of the simulated material in the subsequent hardening stage II matches that of unirradiated materials. This observation provides an a posteriori justification of our assumption that, for understanding slip channel formation, it is sufficient to consider deformation that proceeds in single slip and it thus is typical of hardening stage I.
In contrast, a main motivation for the investigation of slip channels is their probable importance for irradiation embrittlement and in particular grain boundary fracture. This requires an extension of the present model to multiple slip conditions and to 3D polycrystals. Such an extension, based on combining wellestablished dislocation density-based CP models for multi-slip conditions with the present stochastic CP approach, will be pursued in further studies.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.