From Statistical Correlations to Stochasticity and Size Effects in Sub-Micron Crystal Plasticity

Metals in small volumes display a strong dependence on initial conditions, which translates into size effects and stochastic mechanical responses. In the context of crystal plasticity, this amounts to the role of pre-existing dislocation configurations that may emerge due to prior processing. Here, we study a minimal but realistic model of uniaxial compression of sub-micron finite volumes. We show how the statistical correlations of pre-existing dislocation configurations may influence the mechanical response in multi-slip crystal plasticity, in connection to the finite volume size and the initial dislocation density. In addition, spatial dislocation correlations display evidence that plasticity is strongly influenced by the formation of walls composed of bound dislocation dipoles.


Introduction
Crystal plasticity modelling of a macroscopic cylinder typically requires elasto-plastic constitutive laws. Usually, the onset of crystal plasticity is modeled through a smooth, continuous transformation [1,2], even though in the rare absence of pre-existing mobile defects it is a fact that the plasticity transition is discontinuous (see Figure 1). In contrast, during nanopillar compression, mobile defects are suggested to be absent [3,4] and the transition is characterized by discontinuous abrupt event sequences (nanoscale) [3,[5][6][7][8][9][10][11][12][13][14][15]. Naively, one might expect that the averaging of an abrupt nanopillar response would lead to a discontinuous average response at the nominal yield point. However, unconventional size-dependent nonlinear ensemble average behavior emerges during quasi-static nanopillar compression of crystals as size decreases [16,17].
In uniaxial compression of microscopic crystals, discontinuous plastic yielding may be realized by considering a collection of randomly placed dislocation sources (pinned dislocation segments) in an otherwise dislocation-free crystal (see Figure 1). However, even in such an idealistic case, after loading to a finite strain, the unloading process to zero stress will leave a corresponding plastic strain and dislocation structure. Reloading to the flow stress appears quasi-continuous, but the behavior is typically nonlinear and "anelastic" [18][19][20][21][22][23][24], originating in locally irreversible but small deformations that correspond to abrupt jumps of pre-existing dislocations. Experimentally in small volumes, it has been found that uniaxial compression of crystalline nanopillars ranging from ∼100 nm to ∼10 µm is characterized by the absence of mobile dislocation segments ("exhaustion" mechanisms), leading to abrupt events and jerky loading responses [7,[25][26][27][28].
The ensemble average of small-volume abrupt behavior, smooth and nonlinear, resembles macroscale crystal plasticity. For uniaxial compression of cylinders, due to the absence of geometric gradients, it is natural to consider crystal plasticity a a local phenomenon [29]. Thus, it is expected that the ensemble average of nanopillar responses should equal the spatial average response of a macroscopic cylinder. Nevertheless, recent experiments [16,17] displayed strong size dependence for the average mechanical response of copper single crystalline pillars with sizes decreasing from 3 µm to 300 nm, showing increasing curvature during quasi-static loading [1,[30][31][32][33][34]. At which scale does the micropillar statistical ensemble averaged strength and hardening equal the spatially self-averaged ones? (a) An annealed small volume yet with dislocation sources (pinned dislocation segments), loaded to a certain strain, responds abruptly at the elastic-plastic transition (solid:load-control, dashed:displacement-control). If the sample is pre-deformed, reloading (dashed line, reload) shows nonlinearity before reaching the flow stress. (b) The uniaxial compression of a cylinder, at load P, may be thought as a statistical collection of microscale pillars' compression. However, at what scale should each pillar compression be considered for such an averaging to be accurate?
In this paper, we investigate how the statistical ensemble average of plastic, abrupt mechanical response of uniaxially stressed small volumes depends on the system size and pre-existing dislocation microstructure. We perform an explicit but minimal discrete dislocation dynamics model study with one and two active slip systems. Two typical initial dislocation microstructures are utilized: (i) annealed (dislocation free) samples; (ii) "mobile-dislocation-rich" dislocation microstructures created by a prior loading history. We demonstrate that the onset of plasticity and continuous nonlinearity of stress-strain curves is caused by inhomogeneous dislocation microstructures that form under prior multislip loading, composed of dislocation dipoles. We also show that, in this model of uniaxial compression, the very observation of scale free power law avalanche behavior is connected to the emergence of the statistically averaged stress-strain curvature. Based on this model evidence, we conclude that single-slip plasticity may be ensemble averaged by compressed nanopillars with diameters even less than 500 nm. However, multi-slip plasticity may be averaged only by finite volume pillar compression with volumes larger than 2-4 µm.
The paper is organized as follows: Section 2 contains the model description and details of our study; Section 3 is focused on the mechanical response of nanopillars of different sizes and microstructures for multi-slip conditions. In Section 4, we focus on the nonlinearity of statistical ensemble average and its connection to spatial edge dislocation-pair correlations. Avalanche statistics is also discussed for different dislocation densities. In Section 5, we discuss our conclusions in the context of the macroscopic constitutive relations derived by the small-volume response ensembles.

Model Description
The uniaxial compression of a nano/micro-pillar is carried out by two-dimensional (2D) discrete dislocation dynamics, where only edge dislocations are considered in one or multiple slip systems. This is an accurate model for thin films [11,12] and it can be considered as a phenomenologically consistent model for uniaxial nanopillar compression [7,35]. The schematic of the uniaxial compression is shown in Figure 2. Using small strain assumptions, plastic deformation is described through the framework developed in [36], where the material's state determination employs strain/stress superposition. Thus, shape asymmetries related to plastic deformation are effectively not considered. Each edge dislocation is treated as a singularity in an infinite space with Young modulus E and Poisson ratio ν. The application of the dislocation analytical solution, which is valid in an infinite space, needs a smooth image field (ˆ) to ensure that actual boundary conditions are satisfied. Hence, the displacements u i , strains ε ij , and stresses σ ij are written as where the (˜) field is the sum of the fields of all N dislocations in their current positions, i.e., The image fieldsˆare obtained by solving a linear elastic boundary value problem using finite elements, with boundary conditions that change according to the dislocation structure and the external load. Slip planes are spaced at 10b, where b is the Burgers vector magnitude of 0.25 nm. We do not consider slip planes that cross the loading boundaries (see Figure 2) to avoid numerical difficulties induced by dislocations hitting the boundaries. Such assumptions will not alter the plasticity mechanism observed in the sample since the effective slip area is 85% of the sample geometry. The crystal is initially stress and mobile-dislocation free. This stands for a well-annealed sample, yet with pinned dislocation segments left that can act either as dislocation sources or as obstacles. A dislocation source mimics the Frank-Read source in two dimensions [36]. Point obstacles are included to account for the effect of blocked slip caused by precipitates and forest dislocations on out-of-plane slip systems that are not explicitly described. Stress caused by the obstacles is not considered in the model. The strength of the obstacles τ obs is taken to be 150 MPa with 20% standard deviation. Obstacles are randomly distributed over the slip planes with a density that is eight times the source density [35], and a dislocation stays pinned until its Peach-Koehler force exceeds the obstacle-dependent value τ obs b.
A dipole will be generated from a source when the resolved shear stress τ at the source location is sufficiently high (satisfying the condition τ > τ nuc ) for a sufficiently long time (t nuc ). The sources are randomly distributed over slip planes at a density ρ nuc (60 µm −2 ), while their strength is selected randomly from a Gaussian distribution with mean valueτ nuc = 50 MPa and standard deviation 10 MPa. Once the source is activated, a dipole is generated and put at a distance L nuc . The initial distance between the two dislocations in the dipole is at which the shear stress of one dislocation acting on the other is balanced by the local shear stress which equals τ nuc . After a dislocation is nucleated, it can either exit the sample through the traction-free surface, annihilate with a dislocation of opposite sign when their mutual distance is less than 6b or become pinned at an obstacle when the dislocation moves to the obstacle site. Glide is governed by the component of the Peach-Koehler force in the slip direction. For the I-th dislocation, this force is given by where n (I) is the slip plane normal and b (I) is the Burgers vector of dislocation I. The Peach-Koehler force (Equation (4)) includes the stress contribution from all other dislocations in the system (sum off elds) and effective stress (ˆ), considering the external loading and correction fields of the superposition method. Dislocations follow over-damped dynamics, therefore they are driven by their Peach-Koehler forces, and the instantaneous velocity of the I-th dislocation is where B is the drag coefficient. In this paper, its value is taken as B = 10 −4 Pa·s, which is representative of aluminum, with E = 70 GPa and ν = 0.33. Simulations are carried out in an incremental manner, with a time step that is 20 smaller than the nucleation time t nuc = 10 ns. At the beginning of every time increment, nucleation, annihilation, pinning at and release from obstacle sites are checked. After updating the dislocation structure, new stress fields in the sample are determined, using the finite element method [36]. The loading mode is set up to be strain rate controlled at 10 4 /s. We primarily focus on multi-slip loading (two active slip systems oriented at ±30 • relative to the loading direction), and compare our results with single-slip systems, with slip orientation 30 • relative to the loading direction.

The Mechanical Response of Finite Small Volumes in Multi-Slip Conditions
First, we investigate the behavior of samples for fixed total deformation strain (1%) in both annealed and mobile-dislocation-rich samples for multi-slip loading conditions. Pre-existing dislocation microstructures in mobile-dislocation-rich samples are altered through the prior deformation of annealed samples at increasing total strain levels (1%, 5%, 10%), as shown in Figure 3 (example of 10% loading history). If uniaxial compression of "annealed" samples (only dislocation sources initially present) is carried out (cf. Figure 3a), then we find a yield strength size effect and stochastic plastic behavior (cf. Figure 3b) that are qualitatively consistent with experimental findings for uniaxial crystal compression of nanopillars [7,28]. When the developed microstructures at certain loading strain are unloaded to zero stress, a stable dislocation structure of a mobile-dislocation-rich sample forms, as shown in Figure 3c. The mobile-dislocation-rich sample is then loaded to 1% strain as shown in Figure 3d.
The direct comparison of mechanical responses for small finite volumes of different sizes (different colors) and microstructures (dashed vs. solid) is shown in Figure 4. The statistical averages of stress-strain curves based on 50 samples are plotted in Figure 4a. Loading of pre-existing dislocation ensembles to 1% total strain leads to nonlinearity, i.e., a smooth and nonlinear response prior to reaching the flow stress. The average curvature drastically decreases as the sample size decreases i.e., a longer continuous transition from elastic to perfect plastic, in contrast to the expected discontinuous yielding of annealed structures (dashed lines). The observed nonlinear behavior is evidently related to the yield strength size effect in small volumes: while the ensemble average of the yield strength increases as w → 0 (see Figure 4b) for either annealed or loaded microstructures, the yield strength distribution (see Figure 4b) becomes drastically wider with system-size for loaded dislocation configurations, in a qualitative agreement with nanopillar compression phenomenology [37]. By comparing Figure 4a,b, one may notice that the yield stress distribution disparity mirrors the system-size dependence of the anelastic (nonlinear) average behavior. The same exponent that controls the yield strength size effect (σ Y ∼ w −α with α 0.65) [35] is the one that controls the nonlinearity of the average stress-strain behavior (not shown). This finding is consistent with recent observations (e.g., see Ref. [16]-Appendix Figure S4d).

Dislocation Pair Correlations and Single-Slip vs. Multi-Slip Loading Conditions
The effects seen in double-slip loading conditions are not generic. The observed nonlinearity depends on the number of slip systems activated, so we compare results produced in single-slip (oriented at −30 • , see Figure 2) and multi-slip loading conditions. The nonlinearity becomes clear when the stress-strain curve is reconstructed by defining σ r = σ − σ f , where σ f is the flow stress prior to unloading. It is seen in Figure 5, where σ r versus plastic strain p , is plotted that the nonlinearity has a dependence on the sample size for double-slip loading. In contrast, single-slip loading shows no clear size dependent nonlinearity for mobile-dislocation-rich samples (see Figure 5a inset). This apparent discrepancy between single-slip and multi-slip loading indicates a possible connection of this size effect to certain spatial features of dislocation structures that are favorably formed under double slip conditions. Spatial features of dislocation structures may be extracted by: (i) the study of pair correlation functions g ss (r) where s, s ∈ {+, −}, as well as (ii) the sign-insensitive correlation function g(r), with r = (x − x ) 2 + (y − y ) 2 for two dislocations located at r = (x, y) and r = (x , y ). Figure 5b shows g(r) for different w, averaged over 20 realizations. A structural peak forms in g(r) at small distances (∼2 nm, with the slip spacing being 2.5 nm), which signifies the formation of dislocation dipoles. The clustered dislocations are not pile ups (at single slip planes) as we confirmed. The scatter of the pair correlation (errorbar shown in Figure 5b) increases with decreasing sample size, consistently with the variability of yield stress shown in Figure 4b. The origin of the pairs can be traced in the dynamical behavior of the model: Dislocations from different slip systems may mutually approach at a very short distance without annihilation, at the intersection of their respective slip planes. There, a stable structure can be formed by dislocations of opposite signs (see Figure 3c for example). The inset shows the behavior of the average g +− (r): pairs of dislocations with opposite signs are clustered at distances smaller than 3 nm with the peak of g +− (r) being higher as w → 0. Dislocation pairs of opposite-signed dislocations may be viewed as a toy model of dislocation junctions [7], even though such analogies should be considered with care. In single slip loading, as shown in Figure 5c, the peak of the pair correlation function appears exactly at the slip plane spacing 2.5 nm, larger than that in the double slip system. For consistency purposes, we also checked analogous results in samples of different aspect ratios, one of them being shown in Figure 5d for w = 1 µm, and no clear difference is found, thus we conclude that α = 4 is adequate for the purposes of this study.
The very formation of bound dislocation dipoles may not necessarily imply any size dependence of the nonlinear mechanical response [3,6,7,38]. However, the origin of the correlated size-dependent response is indeed tracked down to the stress-field imposed by these inter-slip dislocation pairs. Namely, a single edge dislocation displays a long-range resolved shear stress that has stability lines at 45 • angle with respect to the slip system angle, and an opposite-signed dislocation can combine to form a bound pair at a nearby slip plane. The inter-slip bound dislocation pairs, discussed in this work, apply equally long-range dislocation stress, as the one originating in a single dislocation. The dislocation pair can be regarded as a super-dislocation where the resolved shear stress along the −30 • slip system is plotted in Figure 6a. There are multiple stability lines that can lead to the kinetic formation of stable but weak pairs, ultimately leading to a size-dependent correlated response. For each such super-dislocation, the shear stress sign-changing locations are shown with green lines; along such lines, it is probable to stochastically form a wall of such super-dislocation dipoles. Naturally, the formation of the identified dipoles and the associated patterns should become more probable as the dislocation density increases for the same system size (or as the sample size increases for the same dislocation density). For this purpose, we investigate the effect of different dislocation densities (for w = 2 µm) through creating dislocation ensembles by unloading at different strains (1%, 5% and 10%). We consider dislocation densities (10 14 /m 2 ) that are 2.74, 7.78, 11.1 (see Figure 6c). It is seen that larger initial dislocation density leads to a more pronounced nonlinearity. The g +− (r) is shown in the inset, which signifies that the smaller dislocation density has a higher peak. Our model is benchmarked with experimental data [35], i.e., and it predicts a realistic strengthening size effect and dislocation avalanche statistics in FCC crystals even though at a much higher strain rate than experiments. Assuming a sole dislocation density effect on the strength, we may estimate that the pre-existing dislocation density of the samples in Ref. [16,17] to be 10 13 /m 2 .
The evolution of the average inelasticity may be tracked through the statistics of abrupt events that caused it. The event size statistics is shown in Figure 6d. Event S is the normalized stress drop defined as ∑ i ∈ eventsteps δσ i /σ max where δσ i is the stress drop and σ max is the maximum stress in single realization. It can be clearly seen that the increase of pre-existing dislocation density and inelasticity leads to a power-law behaving ensemble with larger cutoff and decay exponent ∼1.23. For very low pre-existing dislocation density, where crystal plasticity is dominated by dislocation nucleation, one can see that the power-law behavior is almost invisible.

Conclusions
This result is in accordance with a wealth of prior work [7,[39][40][41][42][43][44][45] that have pointed that critical avalanche dynamics requires pre-existing random or "glassy" dislocation microstructures. However, the current work represents a pioneering effort to identify the precise origin of such random structures in small scales. The possible distinction of this work is the fact that realistic dislocation microstructure formation, contrary to a purely random dislocation microstructure, may lead to clear power-law abrupt event statistics and associated effects.
In summary, we identified and studied a nonlinearity in the stress-strain initial-condition ensemble average response during uniaxial compression of small finite volumes. This nonlinear effect is an outcome of small finite-volume avalanche responses [7] and its presence may challenge any possible correspondence between large-scale mechanical response and ensemble averages of small finite volumes (see Figure 1). We find that such correspondence is plausible and sensible for single-slip loading conditions and sample widths down to 500 nm, but not for multi-slip loading conditions with sample widths up to 2 µm. We track the very origin of this effect in the structural features of the emerging dislocation structures and the formation of bound dislocation dipoles.
This dipole formation resembles dislocation junction formation in more detailed models of dislocation dynamics [7]. We may consider a typical macroscopic phenomenological power law strain hardening relation to model this effect in continuum plasticity modeling [1], by stating that the post-yield stress, σ = K n with K the strength coefficient and n the hardening exponent. We find that n (which is defined as logσ log ) is a function of the pre-existing dislocation density ρ, leading to the constitutive relation n = 1 − (644 − 35.35ρ * ) , where ρ * = ρ/ρ 0 with ρ 0 being 10 14 /m 2 . Thus, our explicit discrete dislocation model study of uniaxial compression in small finite volumes demonstrates that ubiquitous abrupt plastic events result into a dislocation density dependent nonlinear dependence.
Together with strength size effects, the identified nonlinearity challenges any attempts for "ensemble averaging" of small-volume responses into forming a representative volume element average. The density dependence can be traced to the pattern formation of microscopic dislocation dipoles, which are not easily formed in single-slip loading conditions. In this way, multi-slip loading conditions are possibly key components to unveiling the role of critical, power-law abrupt events' for phenomenological crystal plasticity.
Author Contributions: H.S. and S.P. designed the study; H.S. carried out simulations. H.S. and S.P. analyzed the data and wrote the manuscript. All authors reviewed the final manuscript.
Funding: This research was funded by the National Science Foundation under award number 1709568.