Microplasticity and yielding in crystals with heterogeneous dislocation distribution

In this study, we use discrete dislocation dynamics (DDD) simulation to investigate the effect of heterogeneous dislocation density on the transition between quasi-elastic deformation and plastic flow in face-centered cubic single crystals. By analyzing the stress-strain curves of samples with an initial, axial dislocation density gradient, we arrive at the following conclusions: (i) in the regime of quasi-elastic deformation before the onset of plastic flow, the effective elastic modulus of the simulated samples falls significantly below the value for a dislocation-free crystal. This modulus reduction increases with decreasing dislocation density gradient: crystals with homogeneous dislocation distribution are thus weakest in the quasi-elastic regime; (ii) the transition towards plastic flow occurs first in regions of reduced dislocation density. Therefore, the overall yield stress decreases with increasing dislocation density gradient; (iii) crystals with dislocation density gradient exhibit a more pronounced hardening stage during which stress is re-distributed onto stronger regions with higher dislocation density until the sample flows at a constant flow stress that is approximately independent on dislocation density gradient. We interpret these findings in terms of a continuum dislocation dynamics inspired model of dislocation density evolution that accounts for inversive dislocation motions. The transition between quasi-elastic and plastic deformation is interpreted as a transition from inversive to non-inversive dislocation motion, and the initial differences in elastic modulus are related to a density dependent polarizability of the dislocation system. The subsequent plastic flow behavior is analyzed in terms of a modified version of Mughrabi's composite model.


Introduction
Spatially graded microstructures can be observed in many biological materials, such as bones, bamboo, shells, etc., where the microstructure change gradually. These graded structures are the result of natural selection and natural evolution, and almost all of them have some special excellent properties [1]. Therefore, inspired by the gradient biological material, metallic materials with graded microstructures were also fabricated and have shown excellent mechanical properties, fatigue properties, friction and wear behavior [2,3]. Gradients of microstructural properties may exist on the levels of grain size distribution, twin density, dislocation density, solute or precipitate density, or combinations thereof [4]. There are many studies of the effects of grain size gradients on plasticity [5][6][7], but there are few studies of the effects of dislocation density gradients despite the fact that such gradients are ubiquitous. Because of the loss of dislocations at surfaces, the existence of a dislocation density gradient in the near-surface region and a concomitant dependency of hardening on distance-to-surface can be demonstrated even in uni-axial deformation of macroscopic samples [8,9]. Even more pronounced are gradient effects in samples where an intrinsically heterogeneous deformation mode imposes a gradient of strain and hence of strain hardening, see e.g. an investigation of dislocation density gradients of pure nickel under torsion, where deformation and dislocation density gradients were found to be accompanied by a decrease in Vickers hardness from surface to center [10,11]. 3 The purpose of our study is to understand the effect of dislocation density gradients on the small-strain deformation of a face-centered cubic (fcc) metals by using the discrete dislocation dynamics simulation method to investigate the deformation behavior in the microplastic regime and during yielding. The paper is organized as follows: We first give a concise presentation of the simulation method and present the results obtained by DDD simulation, which we discuss in terms of characteristic differences of dislocation behavior in the quasi-elastic/microplastic and the plastic deformation regime of microcrystals with and without dislocation density gradient. We then move to a more theoretical description, where we formulate a dislocation density based framework to interpret the characteristic differences between dislocation behavior in graded and non-graded samples, and between quasi-elastic and plastic deformation regimes. We conclude with a brief discussion that puts our results into the general context of recent interpretations of the elastic-plastic transition.

Simulation Set-up
We use discrete dislocation dynamics (DDD) simulation to describe plastic deformation in terms of dislocation motion and dislocation interactions (see e.g. [12][13][14][15]) and apply this method to investigate the effect of dislocation density gradients. Specifically, we use the Parallel Dislocation Simulator (ParaDis) code which uses an adaptive nodal description of the dislocation lines to trace the evolution of the dislocation system [16]. For a detailed description of the code, the reader is referred to the literature [12]. Originally developed by Lawrence Livermore National Laboratory for large-scale hardening simulations [12], ParaDis has become a popular tool for DDD simulation that has been adapted to a wide range of problems including dislocation processes at surfaces [17], size dependent plasticity of small samples [18], and interactions of dislocations with twin boundaries [19]. 4

Initial and boundary conditions
Our initial dislocation configuration consists of randomly located Frank-Read Sources, i.e., initially straight segments with pinned endpoints and a source length l FR = 400 nm. We create a total of n FR = 600 sources in the volume V of each simulated sample, which amounts a total dislocation density of 12 2 76.5 10 m The sources are equally distributed over the 12 possible slip systems. The orientations of the pinned segments within the respective slip planes are equi-distributed over the unit circle, and the segment centerpoints are located at random positions, under the provision that the segment does not intersect the sample surface.
Sources with pinned endpoints represent artificial configurations and it is necessary to ensure that the artificial pinning points do not control dislocation motion and dislocation multiplication. To avoid artefacts caused by the initial pinning points, it is necessary to ensure that the dislocation source length is significantly larger than the dislocation spacing: In that case, operation of a source is controlled by interactions with other dislocations and by the requirement that the emerging loop can move through the dislocation forest. As discussed in [20], this can be ensured by taking the source length to be at least 2.5-3 times the dislocation spacing, i.e., . This condition is in our simulations always fulfilled. It ensures that the flow stress τ f is not mainly controlled by the Orowan stress for source activation, but by the standard Taylor stress representative of forest hardening, where is the Taylor factor which depends on the distribution of dislocations over the different slip systems and is the shear modulus of the material.
For the present study about the effect of dislocation density gradients on the initial deformation behavior, the samples are divided along their long axis into three layers of width 4000 b and volume V/3. Within each layer, a given number n i of dislocation sources are located, hence, the dislocation density per layer is ( Dislocation densities and dislocation spacings in each layer are summarized in Table 1 together with the overall dislocation density gradients of the three sample types.
Realizations of initial dislocation configurations for the three sample types are illustrated in Fig. 1.
Regarding boundary conditions, we assume that dislocations can freely leave through the sample surfaces. This boundary condition as implemented in standard ParaDis (i.e., without surface corrections to the dislocation stress fields) poses certain problems when compared with simulation schemes that fully account for the elastic boundary value problem [21]. While the stress state for uni-axial deformation is correctly represented on average, the same is not true for the stress fluctuations caused by the intrinsic heterogeneity of dislocation slip. In particular, image forces acting on near-surface dislocations are not accounted for. This implies that, in the near-surface region where image forces prevail over the mutual interactions of dislocations, our simulations do not accurately represent the actual dynamics. This region can be 7 estimated to consist of a layer with a width of about μm from the specimen surfaces. As a consequence, the results of our simulations, while representing a qualitatively correct picture of the effects caused by surfaces and dislocation density gradients on plastic deformation behavior, cannot be considered fully accurate in quantitative terms. However, for the case of a homogeneous dislocation distribution, the quantitative differences in comparison with comparable simulations that account for the correct surface boundary conditions [22] are minor.
For the purposes of gaining a qualitative understanding of the phenomena, which is what we are aiming for, the present simplified approach is therefore fully adequate.

Simulation results
In the discrete dislocation dynamics, the initial configuration has a significant influence on the simulation results. Due to the different initial configurations, the stress-strain curves and the dislocation density evolution are different in the presence of a dislocation density gradient. Fig. 2 shows the stress-strain curves for the sample without dislocation density gradient and for the two samples containing a heterogeneous dislocation distribution, one with a small and the second with a larger dislocation density gradient as specified in Table 2. From the figure, we can see that the stress-strain curves exhibit two distinct regimes. In the regime of small strains, we observe an approximately linear increase of stress with total strain. However, the behavior in this stage is not simply elastic. This can be seen from the fact that the slope of the stress strain curves (shown in the inset of Fig. 2) falls significantly below the theoretical value for an isotropic material with ν = 0.324 and μ = 54.6 GPa, for which we expect an elastic slope of E = 2μ(1-ν) = 144.4 GPa. The actual values of the quasi-elastic slope (henceforth: effective elastic modulus E eff ), however, are close to 120 GPa. Moreover, the effective elastic modulus is found to depend systematically on the dislocation microstructure. From the initial slope of the stress-strain curve, we find an effective elastic modulus, which for the homogeneous dislocation arrangement amounts to 118.8 GPa, for the sample with weak dislocation density gradient to 120.4 8 GPa, and for the sample with high dislocation density gradient to 123.5 GPa. We can thus conclude that, in the initial quasi-elastic regime of the stress-strain curve, the presence of a dislocation density gradient leads to a stiffer response of the material.  The opposite behavior is observed during the initial stage of plastic flow. The transition from quasi-elastic to plastic behavior occurs earliest in the material with high dislocation density gradient, later in the sample with small density gradient, and even later in the sample with homogeneous dislocation density distribution. 9 Accordingly, during the initial stage of plastic flow, the presence of a dislocation density gradient makes the material softer. On the other hand, the initial hardening slope increases with increasing dislocation gradient. As a consequence, the flow stress differences between samples with graded and with homogeneous dislocation microstructure decrease with increasing plastic strain, and above a total strain of about 0.35% (plastic strain of about 0.25%) the flow stress of all samples saturates at a level of about 170 MPa that does not depend on the initial microstructure.
The onset of plastic flow is not associated with any marked increase in the total dislocation density, which, as shown in Fig. 3, changes by less than 10% during the entire simulation. The same is true if we look at the three layers individually (inset in This finding is remarkable: The simplest explanation for the observed deformation behavior would be that, in line with Eq. (2), deformation starts in the 10 layer where dislocation density and flow stress are lowest. The ensuing local plastic flow then would lead to an increase in dislocation density that, through isotropic hardening, compensates for the initial flow stress differences such that, ultimately, the sample deforms homogeneously and the flow stress reached no longer depends on initial conditions. This explanation is wrong, as the simulations neither indicate a significant increase of dislocation density with plastic strain, nor a reduction of the dislocation density gradient. We will discuss, in Section 4, the reasons for this behavior and propose a better explanation for the observed flow characteristics and dislocation density evolution.

Analysis and dislocation density-based modelling 4.1 Quasi-elastic deformation regime
To understand the observed behavior regarding the influence of dislocation density gradients on plastic flow, we first consider the initial regime of quasi-elastic deformation. In this regime, the total strain consists of an elastic contribution, plus a strain that is due to displacement of dislocations out of their configurations of minimum energy. Since this displacement is reverted upon unloading, we cannot call this a plastic strain. On the other hand, because the corresponding motions are associated with energy dissipation, we can also not call it reversible in the thermodynamic sense. We therefore adopt a proposal of Zaiser and Seeger [23] and use the neologism 'inversive' for this type of geometrically reversible, but thermodynamically irreversible dislocation motion. Accordingly, the strain caused by inversive dislocation motions is denoted as the inversive strain inv  , and the transition from the initial quasi-elastic to the subsequent plastic deformation regime is understood as a transition from inversive to non-inversive dislocation behavior.
To model the behavior of the dislocation system in the inversive regime, we first consider the paradigmatic case of widely spaced Frank-Read sources that act independently. We use a line tension approximation where we envisage the pinned source segments as elastic lines of line tension where K depends on line direction. For Cu, DeWit and Koehler [24] give an approximately 11 sinusoidal dependency of K on the orientation angle θ with average value 59.3 K  GPa. The shape of the bowed out segment is then an elliptical segment (oblate for screw and prolate for edge orientations). We average over all orientations to obtain an average line tension with and a segment shape that is, on average, circular. The radius of curvature of such a segment is , and in the process of bending the segment sweeps the area The corresponding strain is inversive as long as the source does not pass its critical configuration. Assuming a system of non-interacting sources of volume density , of which a fraction is located on active slip systems with non-zero Schmid factor M, the inversive strain is in the low-stress regime given by We re-write this in terms of the dislocation density, axial stress and Young's modulus where we used that . The total strain in the quasi-elastic regime then follows as the sum of the inversive and elastic strains, from which the effective elastic modulus is deduced as . We thus expect the effective elastic modulus to 12 be reduced due to the presence of inversive dislocation motions.
In quantitative terms, with and we get an effective modulus 115GPa which is close to the values found in our simulation, though slightly too low. We also note that, if we calculate the average inelastic strain from Eq. (7), the resulting effective elastic modulus is independent of the presence or absence of a dislocation density gradient, in disagreement with the simulation results. Thus, the assumption of independently polarized FR sources needs to be modified.
We can arrive at a better description by assuming that the polarizability of the FR sources is reduced by the influence of the surrounding dislocations by a factor P which describes the reduction, due to stresses of surrounding dislocations, of the effective area between the unstressed source configuration and the saddle-point configuration. We thus replace Eq. (6) by the expression We thus find that interactions between dislocation segments of different sources reduce the polarizability as compared to a system of non-interacting sources. For treating the case of inhomogeneous dislocation arrangements, we expand the polarizability around the reference density  : p=0.0425. Inserting this relation into the expression for the effective modulus, we find for an inhomogeneous dislocation arrangement: where the angular brackets denote spatial averages. It is easy to see that, for a heterogeneous dislocation arrangement of a given average density, the mean square density increases with increasing degree of heterogeneity in comparison with a homogeneous system. Accordingly, the inversive strain is reduced and the effective modulus increases. In quantitative terms, we can evaluate the effective elastic modulus of a graded dislocation system as With 0.0425 we obtain the effective modulus values in Table 3, which are in good agreement with the simulation data. We may thus conclude that the observation of a stiffer behavior in samples with dislocation density gradient results from the fact that the inversive polarizability of the dislocation system decreases when the ratio between the dislocation spacing and the spacing of pinning points is decreased.

Plastic deformation regime
The transition towards plastic flow occurs when the externally applied stress 'tilts' the elastic energy landscape to such an extent that dislocations cross barriers in the initial energy landscape, such that upon stress removal they relax into new configurations. In our terminology, this corresponds to a transition from inversive to non-inversive dislocation motions, and from a situation with conserved dislocation density to a situation where both dislocation multiplication and annihilation at the sample surface are present. In order to understand the observed deformation behavior, 14 we thus need to formulate a model describing these processes. To this end, we use the recently formulated approach of continuum dislocation dynamics (CDD) in a version which accounts for dislocation loop generation by FR sources, dislocation line generation by loop expansion, and dislocation losses at the sample surface [26]. We consider the most simple possible CDD version which considers the dislocations on a given slip system in terms of their scalar density ρ and loop density q and note that, because of symmetry, we can treat all four active slip systems in our simulations as symmetry equivalent. Balance equations for ρ and q which account for FR sources have been formulated in [26]. Generalization to include dislocation losses at surfaces is straightforward: Let A denote the cross-sectional area of the specimen volume with an active slip plane. The characteristic rate for a loop to leave the specimen through its surface is given by leading to a loop loss rate and a corresponding dislocation density loss . In our simulations, 1.22 and a simple estimate shows that for these values and our dislocation densities, dislocation loss at surfaces prevails over losses due to mutual annihilation. Combining these relations with the generation rates due to loop emission and loop expansion as formulated in [26] leads to the simple evolution equations The loop emission rate is controlled by the characteristic velocity FR v of dislocations that are crossing the saddle point configuration of a FR source that is interacting with surrounding dislocations. In a bulk material, as discussed in [27], the velocity at the source is indirectly controlled by the velocity v at which dislocations are convected away from the source as these dislocations exert a back stress, which limits source operation. The same is not true in small samples as studied here, since in these samples the emitted dislocations together with their back stress disappear at the surface sink. 15 The dislocation velocity v is controlled by collective dislocation interactions, which we express in terms of a dislocation density dependent 'Taylor stress' . The resulting velocity is given by whereas the source velocity is additionally influenced by the Orowan stress arising from the need to bend the segment between the source pinning points. Assuming that both effects are additive, we have We can now use this framework to look at the steady state that is reached after some straining when dislocation emission and surface losses mutually balance each other. In this case we can look at steady state solutions of Eq. (12). This gives us the following relations connecting, in the quasi-steady state, the plastic strain rate with the initial FR density, the density of gliding dislocation loops, and the velocities of FR source activation and subsequent glide: From this we can infer the FR source velocity in our simulation. We obtain for the which is much less than the dislocation density associated with the source segments.
We can thus conclude, in agreement with the simulation results, that the increase of dislocation density due to the emission of loops is expected to be small, and that the dislocation density remains essentially at its initial level even in the plastic regime.
The reason why the dislocation density does not appreciably increase after the onset 16 of plastic flow, and thus that the initial dislocation density gradient does not disappear, lies in the fact that emitted dislocation loops readily disappear at surfaces and therefore source activation does not greatly increase the dislocation density. What then explains the observation that initial flow stress differences due to a dislocation density gradient disappear in the course of deformation?
To answer this question we note that the flow stress is controlled by the stress needed to activate sources, whereas the rate-dependent contribution to the flow stress is small. We can therefore, for a region that is plastically activated, express the local flow stress by the approximate rate-independent relation which for a homogeneous dislocation distribution gives us, with the parameters in Table 1, an axial flow stress of 175 MPa, in good agreement with the simulation data shown in Fig. 2.
To describe the situation of a gradient dependent dislocation distribution with average gradient , we note that in this case the local flow stress is given where is the midpoint along the specimen axis. To evaluate the corresponding global flow stress, we resort to a variant of Mughrabi's composite model of heterogeneous dislocation density distributions [27]. We assume that the stress is composed of an external stress, which for the present geometry is constant throughout the sample and given by , where the angular brackets denote averages over the sample length, and an internal stress which counter-balances plastic strain heterogeneities. In the original composite model, the Eshelby-like factor η E is set to unity. This is equivalent to making an iso-strain assumption according to which the total strain ε is constant throughout the sample-an assumption which is warranted in bulk deformation but clearly incorrect in axial deformation of micropillar samples as considered here. We therefore relax this assumption by allowing for incomplete elastic compensation of plastic strain heterogeneities, setting η E <1 to correctly reproduce the observed post-yield stress transients. In the plastic regime, we then require the following inequalities to hold: where L p denotes the plastically deforming part of the sample, L e is the part of the sample that undergoes only elastic (or inversive) deformation, and the stress within the plastically inactive region is defined by setting . These relations allow us to evaluate, the mean (external) stress as pl ext el ff  (4) and (6)) where now and the swept area is evaluated according to where now if we are in the quasi-elastic regime, and u=1 in the plastic regime. Evaluating the ensuing local anelastic strains, integrating them over the specimen length, and adding the result to the elastic and plastic strain contributions gives the stress-strain curves shown in Fig. 4 (right). The curves are now in good agreement with the simulation findings. In particular, we observe even for a homogeneous dislocation arrangement a gradual transition from elastic to plastic behavior, which is associated with a divergence of the inversive polarizability: The stress derivative of the inversive strain, Eqs. (21) and (22), diverges at the critical stress where the system passes its critical saddle point configuration. In the stress-strain curves this implies that the stress approaches the flow stress with a horizontal tangent, but not with a discontinuity of slope.

Discussion and Conclusions
We have studied the transition from quasi-elastic behavior to plastic flow in small samples containing either a homogeneous distribution of dislocations or a heterogeneous distribution with a dislocation density gradient. Dislocation motions and dislocation associated strain are already present in the very first stage of deformation, where they are associated with kinematically reversible ('inversive') motions of dislocation segments. In this regime, samples with heterogeneous dislocation distribution exhibit a higher slope of the stress-strain curve. We could interpret this phenomenon in terms of an effective polarizability of the dislocation system that is a decreasing function of dislocation density.
At a critical stress (the flow stress), the polarizability of a homogeneous dislocation system diverges as the system passes over a saddle point and reaches a state of sustained plastic flow. Similar transitions from inversive to irreversibly flowing behavior have been discussed in a wide range of physical contexts, ranging from low-strain deformation of amorphous solids [28] over colloidal suspensions [29], and motion of dislocations under oscillatory stress [30] to vortices in superconductors under oscillating fields [31]. Here we have given for the first time expressions that relate the inversive strain and its diverging susceptibility to basic parameters of the dislocation system in a deforming crystal such as source length and dislocation 20 density.
In the plastically flowing regime we have found that the behavior of In conclusion, we note that the simulation and model predictions compiled in this paper can be easily checked experimentally. Dislocation density gradients are a ubiquitous feature of the near-surface region of deformed samples [8,9]. Owing to general scaling relations, such gradients are bound to grow as the overall dislocation density increases in the course of strain hardening. It is therefore possible to prepare surfaces with different dislocation density gradients and then use FIB to produce micropillar samples for controlled testing. It will be an interesting task for future studies to use this technique to gain a clearer understanding of the deformation properties of samples with heterogeneous and graded dislocation structures. This is particularly desirable in view of the technological importance of processes such as shot peening which induce large dislocation density gradients in the near-surface region in order to enhance fatigue and wear resistance.