**1. 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–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–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–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–34]. At which scale does the micropillar statistical ensemble averaged strength and hardening equal the spatially self-averaged ones?

**Figure 1.** Schematic of ensemble and spatial averaging in uniaxial compression of cylindrical samples. (**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.

## **2. 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 *ui*, strains *εij*, and stresses *σij* are written as

$$
\mu\_i = \mathfrak{u}\_i + \mathfrak{u}\_{i\prime}, \; \varepsilon\_{i\prime} = \mathfrak{k}\_{i\prime} + \mathfrak{k}\_{i\prime\prime}, \; \sigma\_{i\prime} = \tilde{\sigma}\_{i\prime} + \hat{\sigma}\_{i\prime\prime} \tag{1}
$$

where the (˜) field is the sum of the fields of all *N* dislocations in their current positions, i.e.,

$$\mathfrak{u}\_{i} = \sum\_{J=1}^{N} \mathfrak{u}\_{i}^{(J)}, \; \mathfrak{k}\_{ij} = \sum\_{J=1}^{N} \mathfrak{k}\_{ij}^{(J)}, \; \mathfrak{d}\_{ij} = \sum\_{J=1}^{N} \mathfrak{d}\_{ij}^{(J)}.\tag{2}$$

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.

**Figure 2.** Schematic of the discrete dislocation model in this study. Red dots stand for dislocation sources and blue dots represent dislocation obstacles. Black lines stand for slip planes. Slip planes that cross loading edges are not considered to avoid possible numerical issues caused by dislocations pinned at loading edges. This simplification/assumption does not alter the main mechanism of dislocation interactions. The slip system angle are also indicated in the figure separated by the red dashed line.

Slip planes are spaced at 10*b*, 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<sup>−</sup>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

$$L\_{\rm nuc} = \frac{E}{4\pi(1-\nu^2)} \frac{b}{\tau\_{\rm nuc}},\tag{3}$$

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 6*b* 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

$$f^{(I)} = \mathfrak{n}^{(I)} \cdot \left(\mathfrak{d} + \sum\_{I \neq I} \mathfrak{s}^{(I)}\right) \cdot \mathfrak{b}^{(I)},\tag{4}$$

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 of ˜ fields) 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

$$w^{(I)} = \frac{f^{(I)}}{B},\tag{5}$$

where *<sup>B</sup>* is the drag coefficient. In this paper, its value is taken as *<sup>B</sup>* = <sup>10</sup>−<sup>4</sup> 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 104/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.

#### **3. 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).

**Figure 3.** Annealed *vs*. Mobile-Dislocation-Rich dislocation microstructures and uniaxial compression in small finite volumes. (**a**) Initially annealed dislocation structure: large (red) dots stand for dislocation sources, small (blue) dots stand for obstacles, and two slip systems are used. The pillar has aspect ratio *h*/*w* = 4; (**b**) Examples of stress–strain curves of loading-unloading process for different sizes *w*, 20 realizations each; for each *w* three of them are shown. (**c**) Dislocation structure after unloading (one representative structure is shown for *w* = 2 μm), the average dislocation densities (1014/m2) for decreasing *w* are 11, 10.6, 8.6. (**d**) Examples of reloading of the pre-existing dislocation microstructure.

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 (*σ<sup>Y</sup>* <sup>∼</sup> *<sup>w</sup>*−*<sup>α</sup>* 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).

**Figure 4.** "Annealed" vs. "Mobile-Dislocation-Rich" Ensemble Averages in Multi-Slip Conditions. (**a**) Average stress–strain curves of different *w* (50 realizations each). The black dashed line indicates the expected elastic behavior due to the material's elastic modulus. Colored dashed lines are the average stress–strain curves up to 1% strain shown in Figure 3b when loading the annealed microstructure. Solid colored lines are the average stress–strain curves for pre-loaded 10% samples. (**b**) Yielding (defined as 0.1% plastic strain) distribution for different *w*. The line type follows (**a**).

#### **4. 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 *σ<sup>r</sup>* = *σ* − *σ<sup>f</sup>* , where *σ<sup>f</sup>* is the flow stress prior to unloading. It is seen in Figure 5, where *σ<sup>r</sup>* versus plastic strain *<sup>p</sup>*, 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.

**Figure 5.** Single-Slip *vs.* Multi-Slip Loading Conditions and Structural Correlations. (**a**) Reconstructed stress–strain curves for different *w* with *σ<sup>r</sup>* = *σ* − *σ<sup>f</sup>* where *σ<sup>f</sup>* is the flow stress prior to unloading. The inset is the reconstructed stress–strain curves when single slip system is used in the modeling. (**b**) Pair correlation *g*(*r*) of dislocation structure obtained by unloading. The inset shows *g*+−(*r*) correlation. The color of each curve corresponds to *w* show in panel (**a**). (**c**) Pair correlation *g*(*r*) of dislocation structure in double and single slip systems. (**d**) Pair correlation *g*(*r*) of dislocation structure in samples of different aspect ratio *α* and fixed *w* = 1.0.

Spatial features of dislocation structures may be extracted by: (i) the study of pair correlation functions *gss*(*r*) where *s*,*s* ∈ {+, −}, as well as (ii) the sign-insensitive correlation function *g*(*r*), with *<sup>r</sup>* = (*<sup>x</sup>* − *<sup>x</sup>*)<sup>2</sup> + (*<sup>y</sup>* − *<sup>y</sup>*)<sup>2</sup> for two dislocations located at **<sup>r</sup>** = (*x*, *<sup>y</sup>*) and **<sup>r</sup>** = (*<sup>x</sup>* , *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.

**Figure 6.** Spatial correlations and the origin of the stress-response dependence on initial dislocation density. (**a**) Resolved shear stress (unit GPa) along −30◦ slip system from a +− dislocation pair that typically form in the system, functioning as super-dislocations. (**b**) Zoom-in of typical dislocation structure formed after unloading for *w* = 2 μm. (**c**) Reconstructed stress *vs.* plastic strain for structures with different dislocation densities. The inset shows *g*+−(*r*) for different dislocation densities. (**d**) Statistics of plastic events during reloading. Plastic event *<sup>S</sup>* is defined as *<sup>S</sup>* = <sup>∑</sup>*<sup>i</sup>* <sup>∈</sup> eventsteps *δσi*/*σ*max. The red line stands for well-annealed system (see Figure 3a left) loading to 1% strain.

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 (1014/m2) 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 1013/m2.

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 <sup>∑</sup>*<sup>i</sup>* <sup>∈</sup> eventsteps *δσi*/*σ*max where *δσ<sup>i</sup>* is the stress drop and *<sup>σ</sup>*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.

#### **5. Conclusions**

This result is in accordance with a wealth of prior work [7,39–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<sup>n</sup>* with *K* the strength coefficient and *n* the hardening exponent. We find that *n* (which is defined as log*<sup>σ</sup>* log ) is a function of the pre-existing dislocation density *ρ*, leading to the constitutive relation *n* = 1 − (644 − 35.35*ρ*∗), where *ρ*<sup>∗</sup> = *ρ*/*ρ*<sup>0</sup> with *ρ*<sup>0</sup> being 1014/m2. 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.

**Acknowledgments:** We would like to thank X. Ni and E. Van der Giessen for inspiring discussions. We also acknowledge the use of the High Performance Computing System (Spruce Knob) of West Virginia University.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.
