Microstructure evolution of compressed micropillars investigated by in situ HR-EBSD analysis and dislocation density simulations

With decreasing system sizes, the mechanical properties and dominant deformation mechanisms of metals change. For larger scales, bulk behavior is observed that is characterized by a preservation and significant increase of dislocation content during deformation whereas at the submicron scale very localized dislocation activity as well as dislocation starvation is observed. In the transition regime it is not clear how the dislocation content is built up. This dislocation storage regime and its underlying physical mechanisms are still an open field of research. In this paper, the microstructure evolution of single crystalline copper micropillars with a $\langle1\,1\,0\rangle$ crystal orientation and varying sizes between $1$ to $10\,\mu\mathrm{m}$ is analysed under compression loading. Experimental in situ HR-EBSD measurements as well as 3d continuum dislocation dynamics simulations are presented. The experimental results provide insights into the material deformation and evolution of dislocation structures during continuous loading. This is complemented by the simulation of the dislocation density evolution considering dislocation dynamics, interactions, and reactions of the individual slip systems providing direct access to these quantities. Results are presented that show, how the plastic deformation of the material takes place and how the different slip systems are involved. A central finding is, that an increasing amount of GND density is stored in the system during loading that is located dominantly on the slip systems that are not mainly responsible for the production of plastic slip. This might be a characteristic feature of the considered size regime that has direct impact on further dislocation network formation and the corresponding contribution to plastic hardening.


Introduction
Plastic deformation of bulk metals is usually a smooth process due to the large dislocation content of the crystal. As the specimen size decreases to or below the micrometer scale the properties of deformation change profoundly in many aspects. First of all, as it was shown by compression experiments of single crystalline micron and sub-micron sized pillars (called microor nanopillars, respectively), deformation gets localized in distinct slip bands giving rise to an inhomogeneous slip surface [1]. Secondly, these distinct slip bands appear in an intermittent fashion in the form of sudden strain bursts that lead to jerky and stochastic stress-strain curves [2]. Finally, the strength exhibits inverse dependence on the specimen size in this regime, a phenomenon called size effect [3].
The reason for these specific properties of submicron-scale deformation is that the number of dislocations in the volume is rather limited. In bulk samples, upon loading, dislocations move, interact and multiply in the crystal and build up complex structures called dislocation patterns which have a fundamental influence on the mechanical response. In the case of nanopillar compression experiments, however, most of the dislocations can escape the sample before being able to multiply or interact with other dislocations. This process has been termed dislocation starvation [4,5] and results in a limited amount of possible dislocation sources where accumulation of plastic strain is possible. In addition, due to the small specimen size, these sources are typically harder to activate compared to their bulk counterparts, giving rise to the so-called source truncation hardening [6]. Due to the depletion of the dislocation content from the sample strain hardening is often not observed at this scale. These processes have been successfully modelled with discrete dislocation dynamics (DDD) simulations that track the motion of individual dislocation lines [7,8] and the size effects were explained in terms of weakest link arguments [9][10][11].
In this paper we investigate the intermediate, so far less understood, size regime between the bulk and the dislocation starved limits. Here a significant dislocation content is preserved or built up in the sample during deformation but several properties of nanoscale plasticity summarized above still prevail. This dislocation storage regime can be observed for micropillars with diameters larger than ∼1 µm. Here we focus on the deformation of face-centered-cubic (fcc) copper single crystals oriented for multiple slip since in such a case dislocation reactions have the strongest influence and are expected to play an important role in dislocation storage. In bulk it is known for decades that under such conditions a cellular pattern develops [12][13][14] that exhibits fractal character [15,16]. Several experimental methods have been applied to assess how this morphology changes at the micron scale and how that influences the plastic response. From post mortem transmission electron microscopy (TEM) measurements conducted on compressed micropillars Zhao et al. concluded that the transition between bulk and dislocation starved regimes takes place between 1 and 5 µm, where the gradual development of a complex cellular structure is seen as the sample diameter increases [17]. According to a box counting analysis the pattern is self-affine with a fractal dimension comparable to that of bulk samples [15]. With increasing diameter they also observed the decrease of strength and strain burst sizes and the increase of strain hardening [17]. These observations are in line with previous investigations performed on Ni micropillars [18]. To determine the type of the dislocations Kiener et al. applied electron backscatter diffraction (EBSD) and conducted DDD simulations to conclude that a significant amount of geometrically necessary dislocation (GND) density builds up during deformation [19]. In situ microdiffraction experiments of Maaß et al. also reports about the development of an inhomogeneous stress state and a large GND density in Cu pillars under multiple slip [20]. In a similar experiment performed on tensile specimens oriented for single slip Kirchlechner et al. concluded that although slip occurred on the primary slip system, GNDs were generated and stored in a large number on the inactive slip systems [21].
To date, no experimental method is available that could fully resolve the internal microstructure in a micropillar in the dislocation storage regime. During the preparation of TEM lamellae the micropillar is destroyed and some dislocations may leave the thin TEM foil before imaging. X-ray measurements, on the other hand, can be performed in situ and the beam may penetrate through the whole sample, but one can measure only average microstructural properties along the beam. Recently a new methodology of high (angular) resolution EBSD (HR-EBSD) has been established that is capable of determining the GND content to high precision on a surface of a single crystal using a cross-correlation based technique on the recorded Kikuchi patterns [22,23]. In addition, individual components of the dislocation density tensor, that characterizes GND density along with the average local Burgers vector direction, can also be obtained which helps in identifying activated slip planes in the crystal [24,25]. This method has been recently applied by Kalácska et al. to investigate the dislocation structures developing in Cu micropillars oriented for multiple slip [26]. Serial focused ion beam (FIB) sectioning was applied and the HR-EBSD measurement was performed subsequently on each surface in order to obtain information in three dimensions (3d) about the microstructure. A well-developed dislocation cell structure was found in deformed pillars but with significantly lower GND density than that of bulk samples which was found to explain the simultaneous observation of strain hardening and size effects. In summary, the experimental investigations suggest the development of a complex dislocation structure in micropillars in the storage regime with a considerable GND content.
The objective of the present paper is to provide further insight on the dislocation microstructure being developed in the dislocation storage regime. We intend to answer the question whether a significant amount of GNDs form upon deformation (although this is not really expected for uniaxial compression) and if so, in which slip systems these dislocations accumulate. To this end, we perform experiments on Cu single crystals oriented for multiple slip to study the formation of the dislocation pattern in situ at different micropillar sizes. We show that, as suggested earlier, with increasing diameter a gradually evolving cellular structure can be observed with a significant locally varying GND content. The Burgers vector analysis suggests that GNDs primarily form on inactive slip systems characterized by low values of the Schmid factor.
In order to explain this somewhat counter-intuitive experimental result we will also perform computational modelling of the microstructure and its influence on the corresponding plastic response. Several approaches are available on various scales to model the evolution of the internal dislocation structure. Molecular dynamics and DDD simulations have been used extensively to understand the specific features of nanoscale plasticity summarized above [7,8,[27][28][29][30][31][32][33][34][35]. However, limitations in computational power prohibits their use at pillar sizes and strain rates investigated in this paper. One, thus, cannot model the evolution of the system on atomic scales or at the level of individual dislocations but needs to use tools on higher scales. In this paper we employ continuum dislocation dynamics (CDD) simulations that consider the evolution of the system in terms of smooth density fields of discrete dislocations. One of the first CDD models was developed by Groma and co-workers in two dimensions (2d) and was derived using statistical physics concepts from the mobility law of individual straight dislocations [36,37]. Although these 2d models are able to capture dislocation patterning of straight edge dislocations in single slip [38,39], by construction they cannot account for, e.g., dislocation reactions, cross slip and multiplication, mechanisms that are expected to play a key role under multiple slip conditions. In 3d, inspired by the pioneering work of Groma, Hochrainer and co-workers formulated a CDD theory that is based on the kinematics of curved dislocation lines and, as such, takes into account dislocation curvature and the spatial evolution of homogenized dislocation microstructures [40][41][42]. To properly describe the dynamics of the system and the interaction between multiple slip systems, a flux-based discontinuous Galerkin scheme has been introduced in [43] and stress terms related to specific dislocation processes have to be introduced in a phenomenological manner. Recently, such formulations were developed to account for cross-slip and glissile dislocation reactions [44] and dislocation sources [45]. The comparison of the results with DDD simulations showed that the proposed scheme yields a physically correct representation of the interplay between dislocation densities on different slip systems. As such, the 3d CDD model has become capable of being successfully applied to realistic deformation problems, like torsion of microwires [46].
The paper is organised as follows. In Section 2 the experimental methods related to micropillar fabrication, subsequent compression, and determination of the internal GND content as well as the basics of the CDD simulations are described in detail. Results obtained by experiments and simulations on the stress-strain response, size dependence and the internal microstructure are presented in Section 3. This is followed in Section 4 by an in-depth discussion of the dislocation mechanisms found to play a key role at the investigated size regime and their influence on the observed deformation behaviour. The paper concludes with a summary on remaining open questions and an outlook to future research directions.

Considered system
In order to get insights into the deformation mechanisms and the underlying microstructure evolution of fcc single-crystals at the micron scale, this paper focuses on the compression of copper micropillars oriented for multiple slip by a 〈1 1 0〉 crystal orientation as shown in Figure 1.
The micropillars have varying sizes in the a = 1-10 µm regime and an aspect ratio of height:width between 2:1 and 3:1 (the latter being shown in Figure 1). In addition to the crystal coordinate system (represented by the Thompson tetrahedron in Figure 1) we introduce an x y z coordinate system attached to the micropillar as shown in Figure 1. Both experimental and simulation results are to be presented in this frame throughout the paper. The slip systems are named according to the Schmid-Boas notation [47] with letters (A, B,C , D) for the slip planes and numbers for the slip directions.
Considering the Schmid factors, four slip systems {B 2, B 4,C 1,C 3} can be identified to play a primary role for a 〈1 1 0〉 crystal orientation, while the other slip systems are expected to be mostly inactive. Dislocations on primary slip systems may be able to interact with each other by Lomer and Hirth reactions additionally to coplanar and self-interactions. Expecting the same amount of plastic slip on all four primary slip systems would result in a plastic distortion tensor containing only positive e x ⊗ e x and negative e y ⊗ e y components. Consequently, the trace of plastic distortion disappears and the volume is preserved. Considering possible interactions of inactive slip systems, i.e. {A2, A3, A6, B 5,C 5, D1, D4, D6}, with the expected primary slip systems {B 2, B 4,C 1,C 3}, the inactive slip systems can be divided into three groups: {B 5,C 5} sharing the same slip planes, {A2, A3, D1, D4} sharing the same Burgers vectors, and {A6, D6} sharing neither the slip planes nor the Burgers vectors.

Experiments
For the in situ experiments a previously heat-treated copper single crystal sample was used [26], where an orientation was set in accordance with the previous section. A sharp perpendicular edge Figure 1. Representation of the real and modelled micropillars as well as the description of the slip systems and their Schmid factors for an ideal uniaxial stress state (σ y y ). Throughout the paper the results will be presented using the x y z coordinate system attached to the micropillar and shown in the sketch. Note that this is different from the crystal coordinate system which is rotated as shown by the given Thompson tetrahedron.
of [110] orientation was prepared by low-energy Ar ion milling at 6.5 kV, 2.4 mA using a Leica EM TIC 3X polisher. FIB milling was then applied to create square shape micropillars close to the edge in lathe milling position [48] using a Tescan Lyra3 scanning electron microscope (SEM). In order to reduce Ga + ion implantation to the surface, FIB voltages and currents were subsequently reduced from 30 kV, 10 nA to 5 kV, 60 pA as the polishing got closer to the final pillar shape.
This way pillars with approximate sides of 1-1.5 µm, 3 µm, 6 µm and 10 µm were created with negligible tapering and height:width ratio of 2:1 to reduce bending upon loading. Furthermore, several pillars were prepared with the height:width ratio of approx. 3:1.
Mechanical tests were performed using an Alemnis in situ nanoindentation platform (Alemnis AG Switzerland) [49]. The in situ frame allowed HR-EBSD mapping while the pillars were kept under load. The targeted maximum deformation for larger pillars (3-10 µm) was 10%, whereas the smallest pillars were deformed up to ∼5-10% with strain rates smaller than 10 −3 1/s. Conductive diamond flat punch tips ranging from 5 µm to 25 µm in diameter were used during the experiments. Mechanical data collected during the compression experiments were analysed by Alemnis AMMDA software. Appropriate corrections were made on the original load-displacement data, including pillar sink-in correction by the Sneddon method [50].
EBSD diffraction patterns were recorded with an Oxford Instruments Symmetry detector (AZtec v4.2 software) and an Edax DigiView camera (OIM Data Collection v7 software). During EBSD mapping an electron beam of 20 kV, 10-12 nA was used. HR-EBSD evaluation was carried out with the BLG Vantage CrossCourt v4 software. Cross-correlation image analysis was performed on 20 regions of interest chosen from all diffraction patterns. The components of the dislocation density tensor (see below) were calculated by a C++ program written by the authors.
In the paper the dislocation density tensor α [51][52][53] is utilized to characterize the GNDs, defined as where dislocations are characterized by their Burgers vector b t and line direction l t for different t types of dislocations. The sum is over all types of dislocations present in the sample, and t denotes the density of dislocations of type t .
As a result of the HR-EBSD measurement one obtains the nine independent components of the elastic distortion β el . Consequently, the full elastic deformation rather than the lattice curvature is used here to calculate the alpha tensor. Information about the measurement accuracy of HR-EBSD compared to EBSD is given in [54]. The values measured at every single point in a 70°s ample tilt typical for HR-EBSD correspond to an average in a surface of around 40×120 nm 2 and the detected backscattered electrons originate from the first ∼50-100 nm region below the surface [55]. Due to the chosen sample geometry, HR-EBSD measurements could only be performed on the [110] surface of the micropillars (that is, the z = 0 plane), and thus components of β el were only available as a function of x and y (using the coordinate system of Figure 1). The α tensor is linked to the elastic distortion β el as where curl describes the postcurl operator [56], i j k denotes the Levi-Civita symbol and Einstein summation rule is assumed. Since the HR-EBSD measurement only yields data in the x y plane and derivation cannot be performed along the z direction, only three components of the α tensor can be determined experimentally, namely [57]: Another procedure to get a lower bound estimation for the GND density can be done by utilizing an optimization method to minimize the total dislocation line energy (L 1 optimization scheme) [58]. Terms that cannot be measured are set to zero, and only pure edge or screw dislocations with the same magnitude of Burger's vector are considered to be present in the sample. These estimations lead to distinguish edge and screw dislocations based on their energies: where ν is the Poisson number. Furthermore, we define α sq = α 2 xz + α 2 y z + α 2 zz which is proportional to the local GND density [26]. In this article, both L 1 and α sq methods are used in order to study GND density distribution and to perform Burgers vector analysis in deformed copper single crystalline micropillars.

Dislocation density based continuum model
The CDD formulation used is based on the framework introduced in [43]. Following [59], we distinguish between a mobile and a network dislocation density. We take into account multiplication processes including glissile reaction and cross-slip according to [44] as well as interaction due to Lomer and collinear reactions. This is complemented by the homogenized dislocation source model introduced in [46]. Although several parts of the CDD formulation are known from the literature mentioned above, the used formulation including the considered stress interaction terms, dislocation reactions and the mobility law is presented in this subsection for a better readability.
The modeling approach incorporates two coupled problems: the elastic (external) problem, calculating the elastic stress field for a given deformation state, and the internal problem, describing the microstructure evolution for a given stress field yielding plastic deformation. The two problems are coupled via the plastic slip.
We describe the elasto-plastic deformation by an additive decomposition of the distortion tensor Du into an elastic part β el and a plastic part β pl . Considering small deformations, the strain ε comprises the symmetric part of the distortion tensor and can be decomposed analogously into an elastic part ε el and a plastic part ε pl . The stress-strain relation describes physical linearity using the Cauchy stress tensor σ and the elasticity tensor C. It holds: (5a-c) The macroscopic equilibrium equation considering body forces f B in the continuum B is complemented by the boundary conditions at the surface ∂B via given displacements u D for the displacements u on Dirichlet boundaries ∂ D B and given traction t N for the applied traction σn on Neumann boundaries ∂ N B. Dirichlet and Neumann boundaries are disjoint (∂ D B ∩∂ N B = ) and cover the surface completely (∂ D B ∪ ∂ N B = ∂B): The plastic distortion is assumed to result solely from the evolution of the dislocation state on each slip system s ∈ {1, . . . , S}, whereby S describes the total number of considered slip systems Here, it is noted that the order in the notation of the plastic distortion tensor is a matter of convention and both d s ⊗ m s or m s ⊗ d s are used in the literature. Since the HR-EBSD method measures the elastic deformation, the elastic part of the distortion tensor is also used for the α tensor calculation in the simulation according to (2).
The dislocation microstructure is described on each slip system by different dislocation densities, specifying the respective dislocation line length per averaging volume, in addition to a curvature density q s . The volume integral of q s can be interpreted as the number of dislocation loops within the averaging volume. Thereby, the total dislocation density Tot s is subdivided additively into the mobile dislocation density M s , which contributes to plastic slip, and the socalled network dislocation density Net s , which is sessile. The mobile dislocation density can again be additively decomposed into a statistically stored dislocation (SSD) density SSD s and a vector of the geometrically necessary dislocation (GND) density κ s , which forms the GND density GND s = κ s . The vector of GND density can be written based on its screw and edge part: The network dislocation density consists of a density of Lomer junctions L s , which contains the line length of the junction between both end-nodes, and a stabilized dislocation density S s , which captures the line length of the stable dislocation network due to their attachment to Lomer junctions. Here, the density of Lomer junctions is shared with the second involved slip system resulting in a prefactor of 0.5 for the calculation of the network dislocation density: Net s = 0.5 L s + S s . Accordingly the total dislocation density can be decomposed in: Tot The internal CDD variables are illustratively shown in Figure 2.
The evolution of the CDD quantities defining the averaged dislocation state within the continuum includes the flux-based kinematic description as well as internal dislocation reactions involving different slip systems. This may result in dislocation multiplication as well as increase, transfer or decrease of dislocations line length within the system. The homogenized, mechanism based dislocation source model according to [46] considers the production of new dislocations loops. Dislocation sources on the individual slip systems can be activated locally if the effective stress on the individual slip system exceeds the respective critical source stress that depends on the current microstructure. Multiplication mechanisms are considered by cross-slip and glissile reactions, that can lead to a complex dislocation mobility between different slip systems whereas Lomer reactions form sessile dislocation junctions stabilizing the dislocation network and therefore typically impede the motion of the involved dislocations. The latter process is reversible and the sessile dislocation network can become mobile due to unzipping of the Lomer junctions. In contrast, the collinear reactions lead to an irreversible decrease of dislocation density due to an annihilation of collinear dislocation lines sections. The evolution equations of the internal CDD variables read: Here, () react is a shortened notation for specific dislocation reactions: () react = () gliss + () Lomer + () coll . Further, any mechanism that results in the formation of a dislocation reaction is indicated by(), whereas the unzipping of Lomer junctions is indicated by() and } covers the decrease of dislocation density on the reacting slip systems due to the mechanisms. For a more detailed consideration of the individual terms, please refer to the literature.
Nevertheless, it should be noted at this point that the link length of the stabilized dislocations due to Lomer junctions L Lomer , defining the strength of these junctions [61,62], is approximated by a Rayleigh distribution in a certain interval. This interval is based on the scattering of the nucleation radii and the travel distances of the dislocations measured within dislocation network of DDD simulations [34]. Furthermore, the averaged link length L Lomer avg is assumed to scale with the averaged dislocation distance and is used as the expected value of the Rayleigh distribution.
The critical stress for unzipping of Lomer junctions τ crit,Lomer s is based on the bow out stress of the dislocation links of Lomer junctions and scales inversely to L Lomer . The critical source stress τ crit s responsible for the formation of new dislocations within the system comes mainly into play for C. R. Physique -Online first, 3rd May 2021 cases of low dislocation densities and therefore considers the critical resolved shear stress τ CRSS , which is taken from the literature. It holds: Using the closure assumptions introduced in [40], the dislocation alignment tensor A s reads Regarding the system boundaries, we assume that the dislocations can leave the continuum unhindered at the surfaces, characterized by the surface normal n, while no new dislocations are allowed to enter the system over the surface. Accordingly, the boundary conditions are described by Robin-type boundary conditions for the mobile dislocation densities, whereas the curvature density vanishes on the inflow boundaries The yield stress accounts for hardening due to the interaction with forest dislocations from other slip systems within the averaging volume. Here, we consider the yield stress term based on [63] and adapt it following [59] including a material interaction matrix a ss > 0 and the shear modulus µ. The effect of elastic anisotropy is averaged out in large dislocation densities and consequently the isotropic shear modulus can be used, see [64]. The variables describes the index of summation over the individual slip systems. Furthermore, the resolved shear stress on a slip system τ s and the back stress τ b s contribute to the effective stress and are given by: The resolved shear stress is determined by the projection of the stress tensor including dislocation eigenstresses on each slip system following [65]. The back stress covers short-range dislocation stress interaction according to [36,66] depending on a material parameter D, or respectively the Poisson's ratio ν, and the shear modulus µ.

Numerical implementation
The CDD formulation has been implemented into a two-scale numerical framework based on [43] using an own customized version of the parallel finite element software M++ [67,68]. Thereby, a finite element approach is used with hexahedral elements and linear ansatz functions for the displacements in the external problem as well as constant ansatz functions and an implicit Runge-Kutta discontinous Galerkin scheme with full upwind flux for the densities in the internal problems. Simplifying but found to be efficient, the same mesh resolution is used for both scales. The time is discretized by an implicit midpoint rule with a fixed time step.

Experiments
Micropillars were deformed in situ to study how the GND density evolution happens at subsequent loading steps (Figure 3). During the experiment, the deformation process was paused while the HR-EBSD mapping was performed on the surface of the pillars. When the process is paused (for approx. 10-25 min during mapping) stress drops can be observed due to stress relaxation in the material and some minor creep of the piezoelectric actuator and the load cell (as seen in Figure 3 σ( ) curve of the 10 µm pillar). Other pillars were deformed from the initial to final stage without these intermittent mapping pauses (Figure 3 σ( ) curve of the 1.5 and 3 µm pillars). Size related hardening effect is also evident here, as not only the yielding point shifts to higher stresses at smaller pillar sizes, but the steepness of the σ( ) in the plastic regime also decreases.
Based on the acquired stress-strain curves, it is evident that in smaller sized pillars larger strain bursts were detected, as the noisiest curve belongs to the 1.5 µm pillar, in contrast to the rather smooth σ( ) curve of the 10 µm pillar. GND density evolution estimated by the L 1 optimization method is quite similar in all samples, showing a significant increase in ∆ GND upon loading. GND density tends to increase faster in smaller pillars, that is most likely an artifact due to the higher EBSD pattern noise level at smaller pillars. In case of reduced volumes, diffraction patterns have increased noise (see Figure 3 right), that inherently increases the obtained values of GND . A series of identically prepared micropillar compression experiments were performed and the resulting σ( ) curves are presented in Appendix A. Note, that the observations made above for the exemplary σ( ) curves apply for all the measured ones. As expected at this size regime, the curves are of stochastic nature, that is, they differ from pillar to pillar due to the randomness of the initial internal microstructure. Interestingly, the highest scatter in the measured stress values do not seem to decrease with increasing size and are still significant at a = 10 µm. In situ HR-EBSD results are shown in Figure 4 for the 10 µm pillar. The corresponding σ( ) curve can be seen in Figure 3 in black. It is worth mentioning that a flat punch tip with a diameter of 15 µm was chosen for the in situ compression in order to avoid shadowing effect on the measured side and maximize the usable surface for HR-EBSD evaluation. As it can be seen in the supplementary video recorded during the compression, a hard tiny particle attached to the punch tip caused a small impression on the top of the pillar, that acted as a source for dislocations to be created in higher numbers in that region. This however was not expected to influence the deformation in the lower parts of the micropillar.
The first map was taken prior to deformation, showing the initial stress state and dislocation density in the pillar. Components of the σ i j stress tensor in the 0% case show no stress concentration in the sample, as expected. Please note the increased values at the top and bottom part of the pillar. These are only artifacts coming from the proximity of the tip and the FIB milling around the pillar. In the top part, the nearby tip acts as excess volume of material that electrons need to pass through before reaching the detector, that causes the diffraction pattern to decrease in quality close to the tip. At the bottom part, remaining of the FIB milled bulk material can act similarly, reducing the number of electrons that can reach the sample and detector surface. Decreasing the image quality can increase the inaccuracy of the cross-correlation based image analysis, hence causing the off-scale values seen at the sample extremities. One diffraction pattern was chosen from this = 0% map as reference (close to the base of the pillar), and all pixels were compared to this point, even at higher strains. This enabled that the σ i j values are related to the "stress-free" state of the material, making the evaluation unified throughout the experiment.
The second map was taken while the sample was loaded elastically. Naturally, the map taken at 0.3% strain shows no major difference from the unloaded state. No notable GND density increase has been detected.
Continuing the compression process, the third map was taken at the onset of plastic deformation, at 2.0% strain. Increased activity in all σ i j components was registered. As seen in the σ y y component, compressive stresses dominate the whole sample, while GND density increases significantly.
At higher strains, slip activity is noted in all stress component maps, while the average GND density reached its maximum value at the highest load of ≈ 9%. Looking at the surface of the pillar, slip bands clearly show multi-directional slip activation. Interestingly, when the GND density distribution is compared to the visible slip traces on the surface, a directional difference is observed. As slip traces on the surface of the pillar correspond to dislocations that already left the system, the GNDs accumulating inside the material do not necessarily align with these slip planes. Indeed, it can be concluded that in some cases GNDs seem to pile up in directions perpendicular to the visible surface slip bands. Figure 5 plots the accessible α i z components of the dislocation density tensor for all studied strains. Similarly to the stress tensor components, the first two stages show no remarkable features in any of the α maps. Increased dislocation activity is only evident at higher strains in the plastic regime, where a rich GND structure is seen to develop with characteristic features in all three tensor components.

Burgers vector analysis
According to the definition of the dislocation density tensor in (1), the α i z components characterize the average Burgers vector of dislocations being parallel to the z axis, that is, perpendicular to the sample surface. In particular, α xz and α y z components in Figure 5 represent edge components of the GNDs, while α zz represents a pure screw component as shown by the complementary sketch in Figure 5.
The vector constructed from the three available dislocation density tensor components B = (α xz , α y z , α zz ) T , is, therefore, parallel to the average Burgers vector of dislocations parallel to the z axis. Since at a given measurement point dislocations on various slip systems may be simultaneously present, one expects B to be a linear combination of different Burgers vectors with weight factors representing the fraction of dislocations in given slip systems. If, however, B happens to be parallel with one of the Burgers vectors of the crystal one can conclude, that GNDs are dominated by dislocations with that specific Burgers vector.
In this spirit, we define the values where b i is the Burgers vector corresponding to slip direction i (see Figure 1) expressed in the x y z coordinate system. The values a i are, thus, equal to cos(ϕ i ), with ϕ i being the angle between B and b i . In Figure 6(a) plots of the a i are seen and values of ≈ ± 1 (≈0) mark regions where B is close to parallel with (perpendicular to) the given Burgers vector. It is evident from Figure 6 that there are extensive regions with |a i | ≈ 1 with various i values. This suggests, that these regions are characterized by the dominance of one type of Burgers vector. In order to better visualize this finding, Figure 6(b) plots at every measurement point  Figure 3 (black line). A sketch of the Burgers vector directions is also included on the right-hand-side, indicating that the sign of the α i z components can be linked to specific types of GNDs. Secondary electron image shows the slip planes on the studied surface.
the Burgers vector index i max that has the highest |a i | value there, and Figure 6(c) shows the |a i max | values. From the latter one can conclude that in ∼90% of the sample surface the angle between B and b i max is less than 27°and in ∼35% less than 18°. So, as expected, B is not exactly parallel with one of the possible Burgers vectors (as it represents an average Burgers vector in the volume element measured by HR-EBSD) but one can clearly find a dominant Burgers vector in most regions of the sample. According to Figure 6  Burgers vector form a granular-like anisotropic structure with preferred directions parallel with the possible slip directions of the B and C slip planes (see Figure 1). We recall, that the active slip planes are B 2, B 4,C 1,C 3 where most of the plastic slip takes place. Interestingly, regions with dominant Burgers vector type 1 and 3 are mostly aligned along the B plane and vice versa, those with Burgers vector type 2 and 4 mostly along the C plane. This suggests that the GNDs are not accumulated on the active slip systems, rather on inactive ones. The fact that the GNDs are yet aligned parallel to an active slip plane is probably because these GNDs were generated by dislocation reactions or cross slip of dislocations moving in the active slip system. In addition, significant regions with Burgers vectors of type 6 are found which is of pure edge character and do not belong to any of the active slip systems.

Size dependence of the GND density network
As already mentioned in the introduction, the GND density network strongly depends on the size of the pillar. In case of the 10 µm pillar, GNDs in the middle of the sample form a cellular dislocation structure (see Figure 5). If the size of the specimen is reduced, the presence of the surface affects even more the collective behaviour of the dislocations. As it is shown in Figure 7, at the diameter of 3 µm the formed GND network after ∼7% deformation looks quite different from the already presented 10 µm specimen. This pillar was only measured post mortem, after the deformation process had finished and the sample was no longer under load. Although some slip planes can be identified, the cellular structure is completely missing. This is a sign that as the pillar size approaches the 1 µm limit reported in the literature, a transition from dislocation accumulation to dislocation starvation takes place. If the surface of the specimen is examined, Figure 7. Dislocation density tensor components of the 3 µm pillar at = 0% and 7.2% deformation. Corresponding σ( ) curve is plotted in Figure 3 (purple line). A sketch of the Burgers vector directions is shown in Figure 5.
slip planes are still clearly observed but GNDs are stored in significantly less amount compared to the 10 µm pillar as evident from the α sq map in Figure 7.
In Figure 8(a) extensive regions with |a i | ≈ 1 can be identified. From the Burgers vector plots in Figure 8(b) one mostly cannot find clear regions that are characterized by the dominance of only one type of Burgers vector, unlike in the larger pillar. At best, some regions with Burgers vector type 6 are observed, but this (similarly to the 10 µm pillar) cannot be linked to any of the active slip systems. Figure 8(c) enforces our conclusion that a dominant Burgers vector cannot clearly be identified in this specimen.

Simulations
The fcc copper micropillars are modelled by the CDD formulation presented in Section 2.3. The elastic constants are given as C 1111 = 168 GPa, C 1122 = 121 GPa and C 2323 = 75 GPa [69,70] leading to a Zener anisotropy of A = 3.22 [69,70]. The isotropic elastic constants needed for the internal stress calculation are determined following [71] as: G = 40.0 GPa, ν = 0.367. The length of the Burgers vector is given as b = 0.254 nm [72]. The material specific drag coefficient is B = 5 × 10 −5 sPa [73]. The coefficients of the material interaction matrix used in the yield stress term are chosen according to [74] as: a self = 0.122, a Hirth = 0.07, a Lomer = 0.122, a gliss = 0.137, a coll = 0.625. The backstress parameter is chosen as derived in [66] as D = 0.256. The constant for the cross-slip probability term is assumed as β = 10 5 according to [75], the activation volume is V act = 300 b 3 [76] and the stress initializing stage-three hardening is set to τ III = 0.028 GPa [73]. The collision rate constants for the dislocation reactions, cp. [59], are set to c gliss = 2c coll = 0.64 and c Lomer = 640. We study specimens with a 〈1 1 0〉 crystal orientation and with a geometric dimension of a ∈ {3, 6, 10} µm according to Figure 1. The system is subjected to a displacement at the top of the system in axial direction (u D | y=3a = −u y e y ) with a rate ofε y y = −4×10 3 1/s up to a maximum loading of ε y y = −10% . The bottom of the system is fixed (u D | y=0 = 0), the circumference surfaces are considered traction free. The critical resolved shear stress of τ CRSS ∈ {101, 65, 47} MPa is given according to [77]. An initial total dislocation density is given as 2.4 × 10 12 1/m 2 equally allocated to the different slip systems and homogeneously distributed over the specimen. The initial dislocation density is assumed to consist purely of dislocation network density, which consist half of stabilized dislocation density and half of Lomer junctions. The considered numerical system resolution is given by a mesh of approximately 12.3 × 10 3 elements leading to an element width of a/16 and a total number of degrees of freedom of about 10 6 . Taking the Courant-Friedrichs-Lewy (CFL) condition into account, a smaller time step have had to be used for the micropillar of a = 3 µm size. Therefore, the maximum applied load was reduced in this case.

Size dependency
To investigate the influence of the micropillar size on the formation of dislocation networks, in the following we vary the micropillar size. The corresponding stress-strain curves and the evolution of dislocation network density are shown in Figure 9. Considering the micropillar sizes of a ∈ {3, 6, 10} µm, the results of smaller pillars show an elastic-plastic transition at higher stresses with a sharper transition and lower hardening at the beginning of plastic yielding. For the chosen 〈1 1 0〉 crystal orientation the elastic-plastic transition occurs at approximately |σ y y | ≈ 235 MPa for a = 3 µm and |σ y y | ≈ 125 MPa for a = 10 µm. Figure 9(right) shows that the change of the amount of network dislocation density, which is pinned within the system, decreases for smaller micropillars. The microstructure evolution within the first two percent nominal strain  leads to an increase of the network dislocation density of 60% for a = 10 µm, of about 30% for a = 6 µm and a slight decrease of about 10% for a = 3 µm.

Microstructure evolution
To gain deeper insights into the dislocation microstructure evolution, the evolution is exemplarily shown for the micropillar of size a = 10 µm. Figure 10(left) shows the plastic slip evolution over the strain considering its allocation to different slip systems. It can be observed that the plastic deformation mainly takes place on the four slip systems showing the highest Schmid factors, i.e. {B 2, B 4,C 1,C 3}, as shown in Figure 10(left). Finally, 95% of the plastic deformation is caused by the activity of these slip systems. Figure 10(right) depicts the total dislocation density during loading as composition of the different dislocation density quantities considered in the formulation. During the deformation, GND density is formed and maintained in the system. GND density represents a fraction of 26% of the total dislocation density with respect to the considered system resolution by the end of loading. The fraction of SSD density of the total dislocation density is 6% and thus lower than the contribution of GNDs. The remaining fraction of 68% of the total dislocation density consists of network dislocation density. The spatial distribution of different dislocation density quantities at exemplarily selected loading steps is shown in Figure 11.
The contribution of the individual slip systems to the dislocation densities is shown in Figure 12. It can be observed that 32% of the total dislocation density is located on the four active slip systems {B 2, B 4,C 1,C 3} for 10% strain. The remaining 68% of the total density are on the eight inactive slip systems. However, differences can be observed in the allocation of density with respect to different types of dislocation density. For example, about half of the network density (48%) is located on the four active slip systems, while the GND density is almost exclusively on the inactive slip systems. A significant amount of GND density (73%) is located on the slip systems Figure 11. Spatial distribution of dislocation densities: total dislocation density, network dislocation density, and GND density shown for different strain states. Caused by the compression loading, the micorpillar shows a characteristic deformation. The resulting components of the dislocation density tensor are shown in Figure 13. Here α xz and α y z have the highest values and show a gradient along the micropillar height (y-direction) and micropillar width (x-direction), respectively. In addition to the three alpha components of the z-plane, that can be determined by the experiment, the simulation results also provide the remaining components of the complete alpha tensor. Taking into account that it is not possible to avoid slight imperfections in experimental studies, we studied a number of possible imperfections that could have an impact on the results. The simulation results of the considered geometric imperfections of the micropillar are given in the Appendix A in Figure 16 showing the resulting alpha tensors.

Discussion
The objective of the present paper is to get further insight on the plastic material behavior and underlying dislocation microstrucuture evolution in the dislocation storage regime of fcc metals. The plastic deformation of the micropillars occurs mainly on the four slip systems with the highest Schmid factors: {B 2, B 4,C 1,C 5}, as was to be expected. This can be observed in the experiment based on the deformation patterns and surface steps of the deformed micropillars (cp. the secondary electron image of the micropillar in Figure 4), taking into account the slip plane orientation (see Figure 1). This is consistent with the simulation results, that enable the determination of the allocation of plastic shear to the different slip systems (see Figure 10). Based on the superposition of the activity on all four primary slip systems (see Section 2.1), the resulting e x ⊗ e x component of the plastic distortion tensor explains the much stronger bulging of the micropillar in x-direction than in z-direction in the simulation.
During loading, an increasing amount of GND density can be detected in both experiment and simulation for all micropillar sizes, see Figures 3 and 10. This is consistent with the findings of Kiener et al. investigating the compression of single crystalline copper micropillars with a high symmetry crystal orientation via in situ experiments and 2d discrete dislocation simulations [19], and the findings of Maaß et al. using in situ Laue and electron backscatter diffraction to examine the compression of single crystalline copper micropillars orientated for double slip [20]. The amount of GND density might be regarded as unintuitive for a homogeneous uniaxial loading of a single crystalline micropillar and raises the question on which slip systems the GND density is formed. And moreover, how the GND density is stabilized in the system.
The investigations presented in this paper show that the GND density is formed on the inactive slip systems for the compressed micropillars. The experimental results, see Figure 4, as well as the simulation results, see Figure 12, show, that the allocation of the GND density is dominated by slip systems that are not mainly responsible for the production of plastic slip. This observation has also reported by Kirchlechner et al. who considered tension tests of single crystalline copper specimen orientated for single slip and presented analyses of in situ scanning electron microscopy and post mortem µLaue diffraction [21]. The simulation results in Figure  12 also show that there is an amount of GND density on active slip systems at the beginning of loading, but it is not preserved in the system during further loading and can leave the system via the surface. In contrast, the fraction of the GND density on the inactive slip systems increases with increasing load and also as part of an increasing fraction of GND density of the total dislocation density within the system. This can be explained by the fact, that the inactive slip systems have a significantly lower Schmid factor, which is zero for an ideally uniaxial external load. Therefore the driving force on the dislocations to leave the system is also significantly lower and the motion might be additionally hindered by existing network structures. Thus, the GND density on the inactive slip systems can easily be stabilized and preserved in the system. However, it should be remarked that the simulation does not resolve the GND and SSD parts in the network separately. Consequently, the network density could also contribute to the GND density in the system, which is not tracked according to the formulation used.
The Burgers vector analysis of the experimental data in Figure 6 identifies different regions and traces of dominant respective Burgers vectors. This yields a hint at the allocation of the experimentally obtained GND density to the individual slip systems. Obtained as an average Burgers vector this represents the distribution of the GND density over the system as shown in Figure 5. Interestingly, the angle of traces (e.g. with Burgers vector 1) does not match the angle of the slip plane of the primary slip system (e.g. slip plane C of slip system C 1) but rather to the other active slip plane (e.g. slip plane B ). From this we conclude that the GND density, or at least parts of it, is located on another inactive slip system. A superposition of active and inactive slip systems (e.g. B 4 and D6) would be conceivable. Glissile reactions between the moving dislocations on the active slip planes (e.g. B 4) and the dislocations on the inactive slip systems (e.g. D6) can lead to new dislocations with the corresponding Burgers vectors (e.g. D1). However, the inactive slip plane (e.g. slip plane D) would run vertically in the cross section.
The resulting dislocation density tensor in the experiment and the simulation shows qualitatively good accordance in its individual components (cp. Figures 5 and 13). The α xz and α y z components show gradients along the micropillar height (y-direction) and micropillar width (xdirection) respectively. Although there can be determined only three (α i z ) of nine components in the experiment, the simulation shows that the most significant components are included. However, it should be noted that in experiments at this length scale imperfections, such as, e.g., variations in the initial microstructure, misalignments, or unintended deformation states, can easily occur and impede or mislead the interpretation of the results. The imperfections can influence the measurements and the material behavior itself and thus lead to a certain scatter in the experimental results, as exemplary shown in Figure 14 in Appendix A for the stress-strain curves for different pillar sizes. An essential advantage of the simulation is, that we can easily modify system and loading characteristics. So even if we don't know whether or not imperfections occur in the experimental results, a systematic investigation of possible imperfections by simulative variations can enable a better interpretation and strengthen confidence in the results. Comparing the experimental results in Figure 5 with the simulation data of the perfect system in Figure  13, we see, that there is qualitatively a good accordance between the experiment and the simulation. However slight imperfections leading e.g. to different gradients in the results can be better understood by taking into account the variations of the results due to imperfections shown in Figure 16. A slight imperfection causing bending in the system might yield a change from the perfectly horizontal transition between positive and negative regions in α xz (see Figure 13) to an inclined transition as shown for the bending case in Figure 16 and maybe slightly present in the experimental results.
Considering the dislocation microstructure in the system, it can be observed that the mobile dislocations partly leave the system via the surfaces or they interact with other dislocations in the micropillar and form a complex network structure. Which of these two processes is more pronounced is related to the size of the micropillar. The transition from dislocation accumulation forming a cellular structure to the dislocation starvation regime can be observed both in the experiment (cp. Figures 5 and 7) and in the simulation (see Figure 9) as the pillar approaches smaller sizes. In the considered size regime, the dislocation network formation is emphasized with increasing micropillar size for small strain states, as shown in Figure 6, and decreases for smaller sizes as shown in Figure 9. This is not distinctively observable any more for larger strains in the simulation. However, the chosen time resolution in the simulation results might affect this observation since the time step could be chosen too large to capture all possible reactions during the ongoing density evolution. This remains to be investigated.
For larger micropillars the probability of network formation increases due to the longer distance a dislocation has to travel to the surface as well as the lower dislocation velocities that can be observed. The formation of the dislocation network structure contributes to the hardening, which can be observed in the stress-strain curves in both the experiments (Figure 3) and the simulations (Figure 9). This is consistent to the findings of Zhao et al. analysing uniaxial compression of cylindrical copper micropillars with a 〈1 1 0〉 crystal orientation via scanning electron microscopy and scanning transmission electron microscopy [17]. However, considering the quantitative differences in the stress values between the experimental and the simulation results, it becomes clear that some mechanisms are not fully covered in the simulation model for the considered systems. So further adaptions of the used formulations, particularly regarding the reaction and interaction mechanisms with respect to the crystal orientation, have to be pursued.

Conclusions
In this paper, we investigated the dislocation microstructure evolution in the dislocation storage regime of fcc metals. Single crystalline copper micropillars with a 〈1 1 0〉 crystal orientation and varying sizes between 1 and 10 µm have been analysed under compression loading. The analysis comprised experimental in situ HR-EBSD measurements as well as 3d continuum dislocation dynamics simulations. The experimental setup provides insights into the material deformation and evolution of dislocation structures during continuous loading. This is complemented by the simulation of the dislocation density evolution considering dislocation dynamics, interactions, and reactions of the individual slip systems and providing direct access to these quantities.
It has been shown, that the plastic deformation of the material is mainly caused by dislocation activities on the four slip systems with the highest Schmid factor. Here, the loss of dislocation density over the surface leads to observable deformation patterns and surface steps. An interesting finding is the increasing amount of GND density in the system during loading. It has been observed for all micropillars considered and is located dominantly on the slip systems that are not mainly responsible for the production of plastic slip. The stabilization of the GND density on these inactive slip systems can be explained by the lower Schmid factor and the associated lower driving force on the dislocations on these systems. However, the presented results indicate that the obstruction for dislocation motion caused by dislocation networks formation may play a significant role for this finding.
In the considered size regime, the dislocation network formation is emphasized with increasing micropillar size for small loading states. This indicates a transition from a dislocation starvation regime to a regime dominated by dislocation accumulation forming a cellular structure. This is accompanied by longer travel distance of dislocations to the surface as well as the lower dislocation velocities that can be observed for larger micropillars enforcing the probability of network formation and yielding a contribution to plastic hardening.
(contract number: NKFIH-K-119561). The simulation work was performed on the computational resource ForHLR II funded by the Ministry of Science, Research and the Arts Baden-Württemberg and German Research Foundation (DFG). KZ and KS acknowledge the financial support for the research group 1650 by DFG under contract number GU367/36-2.

Supplementary data
Supporting information for this article is available on the journal's website under https://doi.org/ 10.5802/crphys.55 or from the author.
Appendix A. Impact of geometric imperfections Figure 14. Engineering stress-strain curves for different sized square pillars (1, 1.5, 3, 6, and 10 µm side length). 6 µm pillars were bent due to misalignment of the flat punch.
Numerous micropillars were compressed experimentally in this study. However, it is well known that to achieve a perfect system and deformation process without imperfections is a quite difficult task. Any deviation from perfect alignment or minor sample preparation differences might have an impact on the results, e.g., the recorded σ( ) curves will differ from each other. Although pillars were prepared sequentially ensuring good similarity, they can't be a perfect replica of one another and, e.g., the initial dislocation density can vary from place to place. The indenter flat punch tip can have a slight misalignment with the pillar's top surface, which will enhance bending upon deformation. For completeness of the investigations shown in this paper, we show all corrected σ( ) curves of the experimental analyses plotted in Figure 14.
In case of the 6 µm pillars the flat punch tip was not well aligned with the top surface of the specimen, hence the change in the steepness of σ( ) at ≈ 0.005. At the early onset of the deformation process, in case of sample-tip misalignment, the crystal is rotated in order to accommodate the mismatch. Afterwards, the sample's σ( ) curve is modified, as more GNDs are introduced into the system through this geometric imperfection, leading to a less pronounced elastic-plastic transition in the stress-strain curve. It needs to be noted that one of the 6 µm pillars was deformed unlike the other in situ tests. When pausing the deformation of this pillar, instead of keeping the displacement constant (like i.e. the 10 µm pillar case) the load was kept constant, hence the difference in the character of the σ( ) curve's shape. In order to avoid any excess slip caused by this continuous punch tip movement, the subsequent in situ experiments were suspended by the "constant displacement" method.
An essential advantage of the simulation is, that we can easily modify system and loading characteristics. So even if we don't know whether or not there occur imperfections in the experimental results, a systematic investigation of possible imperfection by simulative variations can enable a better interpretation of the results. The effects of various imperfections on the spatial distribution of the dislocation density tensor are studied by simulation in this section. Starting from an ideally straight micropillar as a reference system (system a), the following geometric imperfections are considered, as shown in Figure 15. A linear shift of the micropillar in z-direction (system b), a curved bow-out of the micropillar in z-direction (system c), and a misalignment in the crystal orientation around the z-axis (system d). The linear shifted and curved micropillars (system b and c) show the same offset in z-direction at half height. Additionally, an embedded micropillar (system e) is considered to investigate the influence of simplifying assumptions in the simulation compared to the experiment. It should be noted that the sharp notch in the material in the embedded system is a numerical challenge and can lead to inaccuracies.
The effects of the geometric imperfections on the deformed state and the α tensor field are shown in Figure 16. The qualitative distribution of α zx and α z y components is similar in all systems, but the angles and size of the transition regions change. For the linear shifted micropillar (system b), a characteristic structure is additionally observed in the α zz component. Figure 16. Absolute displacements of the deformed micropillar (scaling factor of two) and all nine components of the α tensor mapped to the initial configuration at a loading of ε y y = −5.7%.