Analysis of dislocation microstructure characteristics of surface grains under cyclic loading by discrete dislocation dynamics

The dislocation structure formation under low-amplitude fatigue in fcc metals for multislip loading conditions is investigated using three-dimensional discrete dislocation dynamics. Tools based on graph analysis, a statistical description of stable dislocation arrangements such as dislocation dipoles and prismatic loops are developed and applied. Upon decreasing the loading amplitude one order of magnitude below the persistent slip band threshold, although qualitative microstructural differences are seen, the elementary features of the investigated defects are the same. A critical number of cycles is required to produce sessile Lomer junctions that stabilize the structure and result in dislocation clustering around them. The crystallographic orientation of the crystal with respect to the loading axis results in different patterns strongly linked to sessile junctions, which are analyzed using spatial correlation functions. The increase in irreversible bulk dislocation arrangements results in roughening of the free surface and increase in surface step heights. Furthermore the crystallographic orientation with respect to the free surface is shown to control the dislocation density evolution combined with the macroscopic Schmid factor.


Introduction
Mechanical fatigue due to cyclic loading remains one of the most challenging problems in materials science. For fcc metals, low-and high cycle fatigue have been thoroughly studied experimentally [1][2][3][4][5][6][7][8]. Fatigue failure in fcc metals is found to be often due to surface roughening with increasing number of loading cycles, leading finally to stress concentration and eventually incipient crack initiation. It is well established that dislocations form persistent slip bands (PSBs) causing surface roughness [9][10][11][12][13]. Irreversible dislocation structures are composed primarily of dislocation dipoles and prismatic loops. Earlier in the 1960s several experimental groups investigated those elementary dislocation arrangements and discovered the important role of prismatic dislocation loops in strain hardening [14,15] and identified cross-slip as a necessary mechanism required for their formation [16,17]. Quantification of prismatic loops in terms of shape and type as well as more precise description of the formation mechanism were given by Feltner [18,19]. More recently, unexpected failure of technical applications in the very high cycle fatigue regime (VHCF) for more than 10 7 cycles and below the PSB-threshold caused to readdress the concept of fatigue limit. Ultrasonic- [20] and piezoelectric [21,22] modern fatigue testing facilities with loading frequencies up to 20 kHz were developed to investigate this problem. Similar to low cycle fatigue, it was shown by [23,24] that fatigue failure in polycrystalline copper in the VHCF-regime is due to cyclic strain localization within PSBs. Nonetheless, details and mechanisms of the microstructural behavior at the dislocation level remain unexplored.
Since dislocation structure formation is controlled by geometrical dislocation arrangement due to elastic interactions between dislocations, the meso-scale simulation method discrete dislocation dynamics (DDD) gives detailed insight into the structure formation mechanisms. Furthermore the role of microstructural parameters, e.g. grain shape and orientation, on structure correlations is accessible. Several groups have used DDD under cyclic loading to explore the mechanisms leading to irreversible dislocation structure formation. A detailed analysis of single-and double slip loading conditions were given in [25,26]. It was shown that an initial single Frank-Read source can yield, due to annihilation and cross-slip processes only, a dipolar dislocation microstructure which shows good agreement with experiments. It was found that slip band spacing was inversely proportional to the square of the slip band thickness. The authors stress the role of cross-slip for strain localization and show that double slip conditions lead to more damage than single slip. Based on statistical distribution of dipole heights, the same group [27] derived a model which could reproduce precisely an experimental hysteresis loop. In [28] results using DDD showed the role of source position, character and length on prismatic loop development while molecular dynamics demonstrated the complex transformations that dipoles undergo to result in debris. Size effect of the plastic behavior during torsional fatigue was investigated in [29]. Recently, it was found that dislocations, starting from a single source, and trapped within a stress gradient can lead to irreversible macroscopic plastic deformation under cyclic bending [30] with very little dislocation storage. The influence of grain morphology and crystallographic orientation were investigated in [31]: surface grains with grain boundaries oriented favoring screw dislocation pile-ups were found to show more irreversible dislocation structure formation than the ones that do not. Therefore, the three-dimensional grain geometry is crucial to understand deformation processes and the strong effect of such special grain boundary alignments have to considered, when setting up model systems. Prismatic loop generation and interaction mechanisms were studied in [32]. The authors present two different transport mechanisms of prismatic loops. If the loops are small in size, then they are driven by stress gradients, whereas long dipoles, e.g. in PSBs, move primarily due to a strong attraction of Modelling Simul. Mater. Sci. Eng. 27 (2019) 055004 T El-Achkar and D Weygand mobile screw dislocations. For multislip orientations, the effect of initial dislocation density and crystal size on surface roughness and dislocation cell size were investigated in [33,34]. The cell size and wall thickness showed an increase upon increasing the crystal size, whereas early surface slip localization was observed for rather small crystals. Several measures such as the Hausdorff dimension for surface analysis and information entropy for pattern analysis were also used. The role of double cross-slip in the development and evolution of slip bands was emphasized.
In this paper, the authors first present methods to quantify irreversible dislocation defects. The new analysis approach is based on a graph description of the resulting dislocation network, allowing for a systematic characterization of the structure. Since direct numerical simulation of high number of cycles is not possible due to high computational cost, the analysis of the microstructure is limited to early stages, giving insight to formation processes and arrangements of stable dislocation substructures. The analysis is applied to dislocation structures of various multislip orientations subjected to two different plastic strain amplitudes. The evolution of the number-and type of defects per cycle shows a two stage behavior within the investigated range of cycles.

Discrete dislocation dynamics (DDD)
The current fatigue simulations are performed using the DDD framework detailed in [35][36][37][38]. Only a brief summary is given here. Dislocation lines are discretized using straight segments joined by nodes. Nodal forces required to move dislocations are computed from mutual interaction with other dislocations and the applied boundary condition. Dislocation reactions are treated as discrete topological events based on energy criteria [36,37,39]. Cross-slip is modeled based on a probabilistic approach [36,40].
As model material, aluminum is chosen and assumed to be elastically isotropic (shear modulus μ=27 GPa, Poisson ratioν=0.347, lattice constant a=0.404496 nm and drag coefficient B=10 −4 Pa s). The simulation volume consists of a cuboidal elastic domain of side length L=2×10 4 a, in which a truncated dodecahedron of volume V dode =1.55×10 12 a 3 representing a surface grain, is embedded (figure 1). The initial dislocation structure consists of Frank-Read sources with a total line length l t and is randomly distributed in the surface grain. The initial total dislocation density amounts to ρ t,0 =l t /V dode ≈4.5×10 11 m 2 . The choice of a surface grain is ascribed to its higher irreversibility compared to volume grains [25]. Moreover, a truncated dodecahedron does not satisfy the crystal symmetry, therefore, unlike cuboidal geometry, exhibits less bias towards screw dislocation cross-slip at grain boundaries [31].
Cyclic loading is implemented by prescribing plastic strain amplitude in the loading/ unloading direction, meaning that the total average dislocation bow out is implicitly controlled. The plastic strain amplitudes refer to the surface grain which is embedded in the cuboidal simulation volume. To analyze the effect of prescribed plastic strain, two different plastic strain amplitudes are used: a large amplitude with ε pl =0.01% and a small amplitude with ε pl =0.001%. The two amplitudes correspond to two extreme cases of high-and low irreversibility for the large-and small prescribed amplitudes respectively [31]. These two strain amplitudes are within the range of other fcc metals and loading conditions, the two limits coincide with the PSB-and irreversibility thresholds given for pure copper in [41] for single slip conditions. In all simulations, the prescribed plastic strain amplitude is symmetric about zero. Boundary conditions are shown in figure 1: dislocations can escape the grain only through the free surface, displacement boundary conditions u y =0 are applied on the bottom surface y=0 and e =˙( ) u Lf t y on the top surface y=L. The applied total strain rate is e =  5000 s −1 and f (t) is a sawtooth function at time t. The sign of the loading rate is reversed once the plastic strain amplitude ε pl is reached. The plastic strain is calculated from the displacement fieldsũ of the current dislocation configuration [36]. The remaining degrees of freedom at the surface have traction free boundary conditions. The elastic superposition problem is solved using finite elements [36,42]. The mesh has a resolution of 32 3 elements using 20-node hexahedral elements.

Graph analysis
A graph representation of the dislocation network allows to employ measures, established in graph theory. The dislocation microstructure is considered as a spatial graph G=(V, E) composed of a set of vertices (nodes) V and edges E representing discretization nodes and segments respectively. A node V i is characterized by a Cartesian position x i , whereas an edge E l connecting nodes V i and V j has properties such as line direction t, i.e. a unit vector pointing from node i to node j, and Burgers vector b acting as a weight. A reaction's Burgers vector b react is calculated using Frank's rule [43,44]. Figure 2 shows examples of reactions leading to a resulting Burgers vector = b 0 react : in (a) the end-nodes of a collinear reaction [45] (figure 2(a)) are connected to each other only by means of virtual edges. The latter type of edges are also stored during the simulation and allow to track the total plastic deformation; in (b) the emitted loop of a glissile reaction (figure 2(b)) with Burgers vector b 3 results in = b 0 react , which means that the junction end-nodes are physically connected to one another only by the glissile loop. Therefore, we differentiate between physical connections as those where the edges carry a non-zero Burgers vector and virtual edges belonging to junctions with an effective zero Burgers vector. Generally, we define the degree of a node in the dislocation graph as the number of connections which a node can have. The following types are identified: Figure 1. Initial and boundary conditions for the surface grain: the grain consists of a truncated dodecahedron and is embedded in an elastic medium: the initial dislocation microstructure consists of randomly positioned and oriented Frank-Read sources, Dirichlet boundary conditions are applied onto both y-ends, other surfaces are set traction free [31]. The application of dislocation concepts in fatigue Dislocations in Solids.
• 1 in the case of pinned end-nodes of e.g. a Frank-Read source or surface nodes. Surface nodes are connected to the outside by virtual dislocation segments to allow for a proper boundary value problem, which requires closed dislocation loops. • 2 which are discretization nodes along the dislocation loops. Those nodes are the most frequent ones and can be added or removed easily during rediscretization. • >2 are related to dislocation junctions. Due to the use of virtual segments, this number is e.g. 4 for a Lomer, Hirth or collinear reaction and 6 for a glissile reaction. The end-nodes of a junction are therefore connected multiple times. More precisely the number of mutual connections corresponds to the number of loops involved in the junction. Nodes that are shared by several junctions may have even more connections.
Based on physical edges, the dislocation microstructure is filtered into components. A component is a subgraph χ, where any two vertices are connected to each other by physical edges and are not connected to additional vertices in G [46]. This was implemented by traversing the edges of vertices recursively and assigning the same number to edges belonging to the same component.

Defect identification
2.3.1. Prismatic loops and dislocation debris. The algorithm of prismatic loop detection is based on filtering the dislocation microstructure into closed components χ c . Among χ c , prismatic loops are subsets which contain at least four distinct dislocation loops, found on at least two glide plane families (with distinct plane normals) 2 and are quasi-planar. The evaluation of quasi-planarity requires the calculation of an average normal c ( ) n c , which is determined using principle component analysis [47]. The positions of m vertices äχ c are written into a´m 3 -matrix Q and centered by subtracting the positions from the average position a c ( ) c . A 3×3 covariance matrix = C QQ T is calculated, and its eigenvalues λ 1 , λ 2 and λ 3 (given in ascending order) are computed. These eigenvalues correspond to the extension of the prismatic loop in different dimensions. In the case of a classical prismatic loop, λ 1 is orders of magnitude smaller than the other two eigenvalues meaning that the eigenvector corresponding to λ 1 is the normal vectorn. A rare case of extremely twisted elongated prismatic loops may show λ 1 ≈λ 2 and l l  3 2 meaning that discretization points are spread rather equally in two dimensions. In this case, eigenvalues are evaluated using the correlation matrix = -- . Finally quasi-planarity is satisfied if the heuristic condition a -< ( ) ·x n a 50 , i.e. the distance from a node äχ c with position x to the plane spanned by the average position a and normaln needs to be smaller than 50a.
Since the discrete dislocation framework assumes dislocation glide, prismatic loops are created by a sequence of cross-slip, collinear-or glissile reactions. Prismatic loops have a Burgers vector which is not contained in the effective loop plane. In the simplest case, they can glide parallel to the Burgers vector direction along prisms of different complexity. In the context of fatigue irreversibility, prismatic loops are particularly stable because the total force acting on them under homogeneous stress fields is zero allowing the loop to only twist due to moments along the glide cylinder [43,48].
The simplest type of prismatic loop consists of interstitial or vacancy type and due to volume conservation of dislocation glide. They effectively have to occur in pairs (one interstitial and the other vacancy type). The pairs are linked by virtual edges within the dislocation network (section 2.2). Since there are several possibilities for the formation of prismatic loops or prismatic clusters, the following classification is done (figure 3): prismatic loops of the first order (figure 3(a)) which arise due to cross-slip and annihilation reactions involving only one Burgers vector. Only two glide plane families (two distinct normals), are involved in this construction and at least four dislocation sectors. This type is simplest prismatic loop. Furthermore, prismatic loops of the second order (figure 3(a)) are defined to comprise three different glide planes and comprise also former Lomer-type junctions (colored in yellow) and non-junction dislocations. At first sight, this looks surprising, as a Lomer reaction generally involves two dislocations. Here only the Lomer reaction and only a single dislocation emerge from the respective end-nodes. This transformation of a Lomer reaction to a spiral source with pinned end-nodes was first observed by [38]. The mechanism leading to this structure involves cross-slip resp. glissile reaction of the two dislocations arms initially forming the Lomer reaction, involving the third glide plane compatible with their respective Burgers vectors. This plane is inclined to the line direction of the Lomer junction. Therefore the incorporated Lomer junctions' end-nodes are pinned, evoke extreme stability and serve as extended obstacles or sources of stress gradients. Among the evolved structures are prismatic loops of equiaxed or elongated shape (figure 3(a)) which have been previously reported in [18,19]. The shape characterization is approached here using the ratio where U is the perimeter and A enc the enclosed area of a prismatic loop spanned byn and a. The value spectrum of the ratio is Î [ ] p 0, 1 shape where p shape =1 for equiaxed and p shape <1 for elongated prismatic loops. The maximum loop dimension for rectangular prismatic loops can be calculated using equations (A.1) and (A.2). The other type of extracted components are classified as dislocation debris (figure 3(c)) which are non-planar subsets of χ c and may contain Lomer junctions. Similar to the second order type, they are very stable.

Dislocation dipoles.
Dislocation dipoles are treated more generally compared to previously discussed defects. The condition to detect dislocation dipoles is purely geometrical, hence, dipolar prismatic loops are a subgroup of dislocation dipoles. A pair of dislocations (dislocations I and II) is considered as a dislocation dipole if both dislocations are parallel, at a distance h<h crit and span a minimum necessary length of l>l crit . More specifically, edges i and j of dislocations I and II respectively are classified as parallel if the closest distance vector connecting the two edges forms an angle larger than 80°with both edges. The critical allowed distance between the edges is calculated using = [48], where τ crit is the critical resolved shear stress required to destabilize a dipole. Here, τ crit is assumed to be the macroscopic total resolved shear stress on the system, which is measured from the stress-strain curve. Substitution yields » h a 35 crit for ε pl =0.01% and h crit ≈47a for ε pl =0.001%. However, due to internal stresses this value can be higher and thus a conservative value of h crit =200a is assumed. The final restriction is that the dislocation pair needs to span a minimum length of L min >2h in order to exceed intrinsic discretization lengths. Furthermore dislocations in the vicinity of grain boundaries (at distances < 100a from a grain boundary) are ignored in this algorithm to avoid mistaking dislocation dipoles with grain boundary pile-ups.

Structure analysis
A fatigued dislocation microstructure is well known to show different types of pattern formation depending on the prescribed strain amplitude, crystal orientation or ease to cross-slip [49,50]. To analyze various resulting dislocation structures, dislocation edges are binned according to their line direction onto ⟨ ⟩ 110 , ⟨ ⟩ 112 and ⟨ ⟩ 122 directions with ≈16°tolerance. In the case of pronounced dislocation bundles, the following method to characterize their spatial distribution is used: dislocation dipoles which belong to a pronounced alignment direction are filtered and the intersection points r i of such dislocations and 200 planes 3 with normal = n t plane are calculated, for which the 2d pair correlation function is computed. In equation (2), r ij is the distance between points i and j, r = On the other hand, microstructures which exhibit a cell structure pattern are analyzed using fractal box-counting technique [51]. This measure has been previously used in the DDD context by [34], however not normalized by the grid size.
A further measure to characterize dislocation density around Lomer junctions is computed using å å where N L is the number of Lomer junction segments, l segment the length of the segment and where d ij is the distance between the middle point of a segment i (belonging to a Lomer junction) and any other segment j. Self-correlation among Lomer junctions is deliberately ignored.

Surface analysis
Extremes at the surface in the form of deep displacement steps are well known to cause stress localization leading to crack initiation. Therefore, to analyze the surface topography evolution during cyclic loading, the surface, which is initially flat, is discretized into a two dimensional equidistant grid with spacing D = a 156.25 . The plastic componentũ z of the total displacement field in z-direction u z is calculated analytically at each grid point. A sensitivity measure to characterize the change in surface step heights during the simulation is the maximum value of the norm of the gradient of the displacement field  ||˜|| u z due to the dislocation motion. The displacement field is evaluated on a regular grid. The gradient x z y z is calculated by computing gradients along the x-and y-directions of the grid using second order accurate central differences. Figure 4. Spatial characterization of a dislocation bundle structure using 2d pair correlation function.

Effect of plastic strain amplitude on microstructural measures for [010] loading axis (setup I)
The role of the prescribed plastic strain is studied in the following section. The initial dislocation microstructure is chosen to be the same for both amplitudes. The initial total dislocation density of ρ t,0 =4.5×10 11 m 2 consists of 40 randomly located and oriented Frank-Read sources with initial lengths between 1000a and 5000a. The crystal orientation with respect to the free surface (normal = [ ] n 001 surf ) and tensile axis ([010] direction) is shown in figure 5(a) (Burgers vectors'-and glide system numbering can be taken from figure A1). The simulations reach 160 cycles for the smaller-and 16 cycles for the larger plastic strain amplitude.
3.1.1. Stress-strain behavior. The macroscopic response of the surface grain is the stressstrain curve which is shown in figure 6 using both (a) the total-and (b) plastic strain. When plotting against the total strain, for the smaller amplitude, the enclosed area is hardly visible. The curve reaches a maximum value of ≈0.085% total strain. The larger amplitude shows a hysteresis with well-defined enclosed areas. Here, the maximum value of total strain reached  Due to the sudden change in loading direction the plastic strain is overshooting. With increasing number of cycles this overshoot decreases and the hysteresis shows a more stable shape. Although the difference in amplitudes is one order of magnitude, the enclosed area of the larger amplitude is approximately five times larger than the one of the smaller amplitude.
The average plastic overshoot is e D » 0.001% o,pl resp. 0.005% for the small resp. large plastic strain amplitude. The width of the overshoot zone is ≈0.003% plastic strain for both amplitudes.
3.1.2. Component analysis. The dislocation microstructure evolution shows for both prescribed plastic strain amplitudes that plastic slip is at first mainly accommodated by the longest dislocations, whereas shorter dislocations act as barriers by forming junctions with mobile dislocations. However, significant qualitative differences in the microstructure are observed. Snapshots of the microstructure colored according to components' number are shown in figure 7, where dislocations colored in black belong to the largest component. For the plastic strain amplitude of ε pl =0.01%, the dislocation density plot shows after 14 cycles a very stable large component alongside with numerous small components that are vastly prismatic loops and dislocation debris. The smaller amplitude ε pl =0.001% shows after about 10 times more cycles a low dislocation density and no clear structure formation. The number of components (figure 8) increases almost linearly with number of cycles, but the slope increases with increasing amplitude by a factor of 50. For ε pl =0.01%, the number of components show rapid increase with a slope change after 6 cycles to reach almost 400 components after 16 cycles, whereas a rather slow evolution is observed for ε pl =0.001% leading to 140 components after 160 cycles. Furthermore, no saturation for the number of components in the investigated range is observed. Figure 9 shows the evolution of the fraction of the component density ρ component with respect to the total dislocation density ρ t for the largest three components. For ε pl =0.01% (figure 9(a)), the largest component constitutes from the early stages and until the end of the simulation the majority of the dislocation density and reaches a saturation value of ≈70% of the total dislocation density after 6 cycles. The component density for ε pl =0.001% (figure 9(b)) shows rather comparable evolution for the largest three components in the early stages and a separation for the largest component after 70 cycles. The largest component fluctuates strongly and reaches peak values of up to 40% of the total dislocation density.
3.1.3. Dislocation density evolution. The evolution of the total dislocation density and dislocation character density are shown in figure 10. Dislocation character is binned according 0.9396 for screw character, whereb is the normalized Burgers vectors. The dislocation density follows a similar two-slope regime as also seen in section 3.1.2. For ε pl =0.01%, the dislocation density reaches almost 16 times higher values than the initial density after 16 cycles while only 3 times higher values for ε pl =0.001% after 160 cycles. Compared to the dislocation density reached after one cycle, the final density for ε pl =0.01% is approximately three times higher. For ε pl =0.001% an 1.2 fold increase is observed. The analysis of the dislocation character density shows that the density of screw dislocations is lowest and fluctuates around a constant value for the examined amplitudes. On the other hand, the density of mixed Figure 8. Evolution of the number of components for both amplitudes ε pl =0.01% and ε pl =0.001%. dislocations is highest and shows a trend similar to the total dislocation density. The edge dislocation density evolution is distinguishable from the rest, in early cycles it is almost equal to the screw density but after a critical number of cycles (6and 70 cycles for ε pl =0.01% and ε pl =0.001% respectively) it shows an increase and separation from the screw dislocation density which is significant for ε pl =0.01% but rather small for ε pl =0.001%.
3.1.4. Defect analysis: prismatic loops, dislocation debris and dislocation dipoles. A detailed analysis of the number of prismatic loops (figure 11) shows the same trend with respect to the slope change which was observed previously for the dislocation density as well as for component evolution. At the end of the simulation, the total number of prismatic loops for ε pl =0.01% is (around 320 after 16 cycles) almost four times higher than for ε pl =0.001% (around 80 after 160 cycles) and no saturation is observed for this value. The larger amplitude shows for the early stages that prismatic loops of the first order are highest in number. After a critical number of cycles, the latter become similar in number to second order prismatic loops and dislocation debris. In the case of ε pl =0.001%, the number of prismatic loops of firstand second order as well as dislocation debris evolve quite similar throughout the simulation. The evolution of the average shape of prismatic loops ⟨ ⟩ p shape during cycling is shown in figure 12(a). ⟨ ⟩ p shape reaches a steady state value of 0.7 for the small amplitude after 20 cycles, while this value is smaller for the large amplitude in the early stages and reaches 0.75 at the Figure 10. Evolution of the total dislocation-and character density for (a)ε pl =0.01% and (b)ε pl =0.001%. Figure 11. Evolution of the number of prismatic loops and dislocation debris with number of cycles for (a) ε pl =0.01% and (b) ε pl =0.001%.
end of the simulation. The distribution of p shape (figures 12(b) and (c)) shows that elongated prismatic loops (p shape <1) occur with a much higher probability for both amplitudes. With increasing number of cycles the maximum shifts from p shape =0.7 to p shape =0.9 for both amplitudes.
The distribution of the side length of the maximum loop extension (see section 2.3.1), which reflects the size of the evolved prismatic loops, is shown in figure 13. For both amplitudes, very small prismatic loops of side lengths a a 50 150 are found more frequently Figure 12. (a) Evolution of the mean prismatic shape. Distribution of the shape factor p shape for different cycles for (b) ε pl =0.01% and (c) ε pl =0.001%.  1 with c 2 =11.5a and c 1 =24.8a resp. c 1 =30a for the large amplitude for cycles 7 and 15 respectively. The maximum dipolar height is located at the closest glide plane distance of h=11.8a. For the small amplitude the fitting coefficients are c 2 =11.5a and c 1 =22.4a resp. c 1 =24.3a for cycles 80 and 160 respectively. Again, the height with the highest probability is very close to the smallest glide plane distance. Increasing the number of cycles in both amplitude shows a decrease in the maximum values towards an increase in larger heights.
The mean length of the dipoles ⟨ ⟩ l ( figure 15(a)) converges for both amplitudes towards a value of ⟨ ⟩ l ≈ 250a. The frequency distribution of dipolar lengths (figures 15(b) and (c)) shows maximum values at l=200a for both amplitudes with a significant drop in this value afterwards. However, it is noticeable that long extended dipoles = l a 1000 2700 exist, although with a smaller frequency. 3.1.5. Spatial arrangement for ε pl ¼ 0:01%. As seen in section 3.1.2, the microstructure developed under large plastic strain amplitude (ε pl =0.01%) is complex and consists of one very large component and many very small components, that are mainly dislocation debris and prismatic loops. The latter were studied in section 3.1.4. A spatial analysis of the dislocation arrangement is done by filtering the edges of the dislocation structure along ⟨ ⟩ 110 , ⟨ ⟩ 112 and ⟨ ⟩ 122 directions (details in section 2.4). In figure 16(a), the fraction of the dislocation density in each of these directions is shown. It is observed that after 4 cycles the  To understand, whether this is a key feature of the system or a specific feature of a single simulation, several simulations were conducted for systems of same initial total density but with dislocations of same initial lengths of either 1000a or 7000a. These results shows that highest alignment is either along Further analysis of the developed structure shows regions of high-and low dislocation density. In order to quantify this pattern using the method described in section 2.4, the simulation domain is discretized into boxes with side length of 400a and only dislocations of lengths above 500a are kept. This allows to remove debris that do not belong to the main structure. Afterwards, a 2000a thin TEM-like slice of the microstructure is taken perpendicular to the [100]-direction ( figure 17(a)). Then the box-counting technique is applied to the dotted square region which contains the largest-and second largest cells denoted by the numbers 1 and 2 respectively ( figure 17(b)). The scaled box-counting measure D ( ) N x 2 shows that slope of this measure approaches zero for box sizes of D » x a 900 and Δx≈10a. In the quasi-linear regime, the slope of this measure is m=45071a.  contour plot ofũ z after 14 cycles during unloading ( figure 19(a)) shows a block-like structure of the displacement field with residual extrusions and intrusions which remain even after unloading. A plot along the diagonal direction ( figure 19(b)) magnifies the nature of the surface displacement; the intrusion is followed by an extrusion giving it a step like character. (laboratory frame) as shown in figure 5(b). The final distribution of dislocation density per glide system after 5 cycles is shown in figure 20 for setup I and II. For setup I dislocation density is distributed among three different types of glide systems. Glide systems 2, 6, 7 and 11 show fairly high dislocation density (almost 17% of the total dislocation density) followed by glide systems 1, 5, 8 and 12 which     Figure 22(b) shows the mean pair correlation function (section 2.4) characterizing the structure: dislocations, which are aligned alongb 3 , are filtered out. The shell thickness and maximum radius are set to dr=30a and r max =4100a respectively. At very small distances, the mean pair correlation function shows a very high value followed by a drastic drop ( figure 22(b). Several well-defined peaks are also observed at r≈700a and r≈1800a.

Stress-strain behavior of the surface grain
The stress-plastic strain curve shows a hysteresis loop for both amplitudes indicating that irreversible dislocation motion occurs leading to energy storage in the resulting dislocation structure. Furthermore the plastic strain overshoots the prescribed values. There are two reasons for this behavior: the loading is applied in a sawtooth fashion and the sign of the imposed strain rate is changed once the absolute value exceeds the imposed plastic strain amplitude. This sudden change and the high applied strain rate give rise to a number of dislocations which can glide further as the stress level is still sufficiently high. Note that in VHCF experiments [21,22], very high frequencies in the kHz regime are used to reach the required number of cycles. Therefore, this effect might be present in experiments too, but this is very difficult to assess as the time resolved deformation and loading conditions of single surface grains are not accessible experimentally.
To quantify the magnitude of the overshoot in plastic strain, a corresponding number of dislocations sweeping a cross-section in the middle of the grain is calculated: 4 dislocations (10% additional dislocations) for the larger amplitude and 1 dislocation (2.5% additional dislocations) for the smaller amplitude can produce the corresponding plastic overshoot. The observed overshoot in plastic strain will therefore not change the fundamental behavior.
The effective stress amplitudes of about 55-75 MPa and total strain amplitudes in the range of 0.085%-0.12% reached in the simulations are both within the range of experimentally studied fatigue life curves for coarse grained and ultra-fine grained Al [52].

Dislocation network characteristics from the graph analysis
Dislocation microstructures under fatigue conditions obtained by different plastic strain amplitudes show at the local scale common features, like the formation of dipoles, prismatic loops and the role of Lomer junctions serving as anchor points for dislocation accumulation. This is characterized by a number of statistical measures with similar evolutions but different absolute numbers, e.g. more or less pronounced two stage regime, with an incubation number of cycles followed by a linear increase with number of cycles in the investigated range.
On the overall scale, the developed dislocation microstructure (figure 7) is qualitatively different: for large amplitude, structure formation at early stages is observed, occupying a large fraction of the total grain volume. Up to 70% of the dislocation density belongs to the largest component while the remainder is found in prismatic dislocation loops and dislocation debris ( figure 9(a)). Furthermore the structure is cellular meaning that there are regions of very low dislocation density surrounded by regions of high dislocation density. The main building blocks of the largest component are very stable dislocation dipoles and Lomer junctions. At this stage, even eliminating artificial Frank-Read sources by connecting pinned end-nodes to each other and thus getting rid of the pinning restriction shows almost no difference during further cycling of the dislocation network compared to keeping sources. This indicates extreme stability of the evolved dislocation network and its independence on initial dislocation sources. Applying this elimination procedure at very early stages results in a large number of prismatic loops which remain stable during further cycling and almost no increase in dislocation density occurs. On the other hand, in the small amplitude simulation, the low dislocation density is contained in the increasing number of prismatic loops and debris clusters. The largest component containing a large number of Lomer junctions shows early stage clustering.
The formation of cellular structure for large imposed plastic strain amplitudes and clustering of edge dislocations for small amplitudes is a known phenomenon [53]. Also the evolution of the number of components (figure 8) shows a nonlinear relationship upon decreasing the prescribed amplitude. Indeed, after 10 times more cycles, the number of components for the small amplitude is approx. 2.8 times less than for the large amplitude. This phenomenon is observed also in additional microstructural measures such as the dislocation density (figure 10) and total number of prismatic loops ( figure 11). Both simulation types show that after a critical number of cycles (6and 70 cycles for the large-and small amplitude respectively), there is a slope change in the analyzed measures. The change is very clear in the large amplitude simulation and less pronounced for the small amplitude.

Dislocation dipoles characteristics and role of junctions
Splitting the dislocation density into different characters (figure 10) reveals also a turning point for the evolution of the edge dislocation density. For ε pl =0.01%, after 6 cycles the rapid increase in edge density is reflected in a microstructure of mean spacing r » a 1 1500 t . A detailed analysis of these edge dislocations reveals that the edge density mainly constitutes dislocation dipoles and Lomer junctions. For ε pl =0.001%, the growth in edge density after 70 cycles is rather slow and related to clustering around the largest component. The analysis of prismatic loops and dislocation debris ( figure 11) shows that prismatic loops of second order and dislocation debris are extremely stable and present sites of dislocation clustering.
The very large number of prismatic loops for ε pl =0.01% is reflected in the smoothness in average values and distribution of the prismatic shape factor p shape (figure 12) compared to large oscillations at early stages for the small amplitude due to the small number of prismatic loops. For the large amplitude, the averaged shape factor ⟨ ⟩ p shape increases with number of cycles, while for the small amplitude it converges to ⟨ ⟩ p shape ≈0.7. From the distributions of ⟨ ⟩ p shape (figures 12(b) and (c)), the stability of the mean value for ε pl =0.001% can be seen since the distributions are very similar with the exception of maximum shift from 0.7 to 0.9 corresponding to prismatic loops of aspect ratios 6 and 2.5 respectively and increase of equiaxed loops at the expense of elongated loops. These results in the early stages are in agreement with experimental findings reported in [18] (aspect ratio of 5). In the case of the large amplitude, there is a drastic shift from a Poisson-type distribution to a cumulative-type distribution and a similar shift in maxima from 0.7 to 0.9. The latter is accompanied by a huge increase in frequency for equiaxed prismatic loops (from about 0.001 to about 0.16), leading also to an increase of the mean value. Furthermore it is found that a large number of loops formed at later stages are extremely small having side lengths in the range of a 50 -100a for both amplitudes.
The analysis of dislocation dipoles shows that the average dipole height ⟨ ⟩ h increases for the large amplitude and remains rather constant for the small amplitude ( figure 14). The distribution of dipolar heights h is exponential, indicating also that the smallest possible distance-limited by minimal distance between parallel glide planes-is the most likely and stable one. The distribution is in contrast to the normal distribution seen in [27] but very similar to the distribution of dipolar heights in PSBs presented by [54]. In the latter work, at the lower end of the distribution and prior to the most probable height, the histogram decreases smoothly to zero due to dislocation annihilation mechanisms, not included in the DDD framework. The authors suggested that dislocation annihilation is not occurring at a critical interplanar distance and fully spontaneous but involves thermal activation leading to dislocation climb. In the present study dislocation climb is not considered and therefore the latter mechanism cannot occur. The maximum height in the experimental distribution at the highest investigated temperature (750 K) agrees with the observations of the simulations. The stability and strength of the dipoles can be understood because of h max =11.8a is almost 3 times less than the critical height h crit =35a, for which a dipole would break apart under a total resolved shear stress of τ crit ≈33MPa. For higher number of cycles the maximal strength decreases and therefore favors larger heights. The average dipolar length ⟨ ⟩ l shows similar convergence in the small amplitude case but decreases for the large amplitude to reach 250a ( figure 15). Since elongated prismatic loops are a subset of dislocation dipoles, the distribution of dipolar lengths l shows strong similarities with the maximum loop dimension. In the models presented by [16,17] short dipoles are a consequence of high cross-slip activity, large strains or double slip orientations. The influence of cross-slip was also demonstrated in [55,56], where dipoles in aluminum were shown to be generally shorter than those formed in copper due to the ease of cross-slip in aluminum.
However, cyclic loading with small imposed strain amplitudes are shown to lead to short dipoles, too. Also long dipoles (>2000a) with a smaller frequency are observed for the large amplitude. Such dipoles are chopped into shorter dipoles by the mechanism presented by [57] involving cross-slip. A detailed alignment analysis of the three-dimensional dislocation network (section 3.1.5) shows, that the cellular structure is composed of dislocations aligned in specific directions. As a matter of fact, most of the dislocations are aligned along ⟨ ⟩ 110 directions ( figure 16). Further analysis shows that Lomer junctions, which are aligned along ⟨ ⟩ 110 directions, are very resistant to cycling and many of them survive throughout the entire simulation. The observation of high strength and stability of Lomer junctions was also reported [58]. These sessile junctions are also present in second order prismatic loops, dislocation debris as well as large complex components. The high dislocation density seen along ⟨ ⟩ 110 directions consists, however, of at most 30% Lomer junctions at the end of the simulation ( figure 23). This observation leads to the hypothesis that dislocation density is higher around Lomer junctions than random dislocations, due to the stability of Lomer junctions. Indeed, figure 24(a) shows the excess dislocation density around Lomer junctions compared to regular dislocation (according to equation (3)). An excess in dislocation density of ≈27% is found around Lomer junctions for r=200a with a drop in this measure for larger radii. DDD simulations [59] for tensile tests on pillars along the [001]-direction reported that dislocations align along all ⟨ ⟩ 110 directions in the (111) plane. In our current study, dislocations alignment shows a bias for certain ⟨ ⟩ 110 directions ( figure 9(b)). To understand the latter phenomenon, both final glide system distribution (figure 20(a)) and final distribution of Lomer junction density per junction Burgers vector with respect to the total Lomer junction density ( figure 24(b)) are examined. The dislocation density among glide systems is split into three categories: type I glide systems (2, 6, 7 and 11) with highest dislocation density followed by type II systems (1, 5, 8 and 12) with medium dislocation densities and finally type III systems (3, 4, 9 and 10) which are almost inactive. From figure 24(b), it is clear that Lomer junctions with Burgers vectors b 1 and b 6 are with the highest density (they correspond to directions parallel to b 6 and b 1 respectively) followed by smaller densities in Burgers vectors b 3 and b 4 and finally b 2 and b 5 . b 1 -and b 6 Lomer junctions result from the reaction of systems of type I and II, whereas b 3 and b 4 from type I and III and b 2 and b 5 from type II and III (the combination of glide systems to yield Lomer junctions of certain Burgers vectors is Figure 23. Fraction of non-junction dislocation density in ⟨ ⟩ 110 direction with respect to the total dislocation density in ⟨ ⟩ 110 direction for ε pl =0.01%. given in figure A2). So the underlying dislocation density distribution correlates with the distribution of Lomer junctions. This means that for the given [010] tensile axis orientation, Lomer junctions arise naturally along [101]and [1̅ 01] directions, which are directions that show highest dislocation density as seen in figure 16(b). These directions are perpendicular to the tensile axis, thus experience no resolved shear stress. Taking the above aspects into consideration, this means that a critical number of cycles is required for stable Lomer junctions to form, which accumulated during the simulation along directions correlated with the dislocation density distribution. Eventually, regular dislocations moving in the vicinity of such junctions get trapped into the structure is aligning accordingly. To characterize possible intrinsic scales of the overall dislocation structure, the boxcounting technique (figure 17) is employed allowing to determine numerically the largest cell size of about 900a, which is related to the drop in the measure at the upper value. The linear domain is extended over two orders of magnitude as in [51] and is related to all other smaller cell sizes. At the lower end, the plateau indicates the wall thickness. But since dislocations are discretized into single points, some of which are isolated, then the measure is related to the smallest resolution of such points.

Surface roughening: role of the orientation of the surface normal
Roughening of the free surface can be seen when examining the distribution of surface displacements ( figure 18(a)). The distribution in the early stages shows that the surface is rather flat with a small standard deviation. In the later stages of the simulation, the distribution is wider (standard deviation increases) and the maximum is shifted away from zero. Although the number of vacancy-and interstitial prismatic loops are equal, they result in extrusions and intrusions at the surface. Critical features for stress concentration, such as step heights, are measured in figure 18(b). It is obvious that this measure increases linearly which means that the accumulation of volume irreversibilities affect the surface evolution.
The dislocation density in setup I (figure 20(a)) is distributed not according to Schmid factor distribution. Although all eight systems of types I and II are symmetric with respect to the tensile axis (Schmid factor S=0.41), it turns out that dislocation density is highest when S surf of the considered slip system is high. Vice versa, systems with small ¢ S show smaller dislocation densities. In other words, the free surface introduces an additional degree of freedom: it is observed that slip systems with Burgers vectors that are parallel to the surface normal show pronounced activity. This means that for type I glide systems ¢ = = S S 0.41, type II ¢ = S 0 and S=0.41 and type III ¢ = = S S 0. A rotation of the crystal system about the tensile axis (setup II in figure 5(b)) yields a dislocation density distribution according to Schmid factor distribution ( figure 20(b)). Although this is still a [010] crystal orientation, all theoretically active systems are symmetric with respect to the surface and, therefore, behave alike. The increased activity of slip systems results in higher volume irreversibility which can be seen e.g. in the increased number of Lomer junctions ( figure 21(a)). Also the emerged surface shows higher roughness ( figure 21(c)) since all eight systems contribute equally to surface steps. The structural difference regarding block-and triangular patterns on the surface is purely a geometrical consequence of the intersection of glide planes with the surface.
To understand the physical reason why the orientation of the free surface impacts dislocation activation, consider the plastic strain tensors for dislocations of Burgers vectors b 4  Lomer junctions are a signature of dislocation structures under multislip cyclic loading which lead to irreversible dislocation accumulation and motion. The structure formation can be understood and predicted by knowing initial-and boundary conditions of the setup as follows: whenever the combination of two systems which have the highest Schmid factor results in Lomer junctions, the microstructural evolution will be governed by the sessile junctions. This result is supported by the experimental observation of Lomer locks in [001] copper single crystals during cycling, which show more rapid hardening than any other crystal orientations. The authors [61] report that due to errors in processing and testing, it is experimentally impossible to obtain an exact [001] orientation, therefore a primary slip system would be initially established. Since in the present ideal aligned setups this is not observed, the sensitivity of the microstructure to a 10°misalignment from the [010] axis is studied. The simulation results also show the presence of Lomer junctions, however, with a reduced role and a rather planar microstructure is observed. Another experimental result related to Lomer junctions was the work of [62] which reported Lomer junctions in a [1̅ 12] double slip orientation. Such junctions are the natural result of slip systems with the highest Schmid factor. In previously conducted simulations, high symmetry orientations showed, compared to single slip, that stable structures form at very early stages in multislip conditions, but require more cycles in the case of single slip. For example, the setup from the classical work of [63] with a 10°away from [110] crystal orientation was simulated. To accelerate the simulation, only the glide system with the highest Schmid factor was considered as the initial dislocation microstructure. The results are in agreement with the experimental work, where elongated dipoles are observed along [12̅ 1̅ ] direction (figure A6) in the crystal frame of reference. However, unloading the structure and eliminating Frank-Read sources using the technique explained above results in dipoles exiting the free surface. Performing the same simulation with all glide systems populated shows that presence of Lomer junctions stabilizes a structure even after unloading and removing Frank-Read sources. It also means that Frank-Read sources act as extra points of structure stabilization. A single glide system in a single slip orientation can reproduce the main part of the evolved dipolar microstructure, however, cellular or more complex structures cannot be mimicked, hence, a minimum dislocation density on other glide systems is necessary to allow for stabilizing dislocation junctions.

Conclusion
In this paper, initial fatigue behavior in multislip is investigated for two plastic strain amplitudes (ε pl =0.01, 0.001%) using three-dimensional DDD. A statistical foundation, based on mathematical characterization of stable defects such as prismatic loops, dislocation dipoles and dislocation debris is established relying on techniques from graph theory. The main outcomes are: • The observed defect characteristics agree with previous experimental observations. Additionally, it is found that the nature of the defects e.g. type, size and shape of prismatic loops, dipolar heights and lengths are qualitatively very similar for the two applied plastic strain amplitudes. The main difference is the evolution rate of such defects with a nonlinear transition between both amplitudes. • The obtained structural patterns for [010]and [110] tensile axis orientations are cellular and slip band structure respectively. Both structures are quantitatively analyzed. The slip band structure, which is a natural outcome of the simulation, shows well-defined spatial correlations firstly observed for multislip DDD simulations. The structure formation is explained by extremely stable Lomer junctions, to which other dislocations attach or align by elastic interaction to form bigger clusters. • A strong influence of the crystal orientation with respect to the surface and not only with respect to the global loading direction is observed and explained: the dislocation density evolution is found to be controlled by both the Schmid factor and the orientation of the Burgers vector of the glide system with respect to the free surface. Glide systems are found to be highly active only if both their Schmid factor and the modulus of the scalar product between Burgers vector's direction and surface normal are large. This means that the full crystal orientation is required (with respect to the tensile axis and free surface) to understand damage evolution within surface grains (for example in the work of [64]). These findings are of importance for the analysis of polycrystalline samples. One of the key issues for fatigue damage is the identification of possible grains more susceptible to damage initiation. In a polycrystal material, assuming a surface with one hundred randomly oriented grains, roughly four grains would lie in that surface orientation spectrum which would be more critical for damage accumulation. The result of volume irreversibilities are also visible on the surface. Both the surface roughness and the step height increase with further cyclic loading. A2. Glide system numbering and combinations to form Lomer junctions Figure A1 shows the numbering scheme for glide systems. The combination of glide systems to form specific Lomer junctions is given in figure A2.

A3. Evolution of the stress peak values
The evolution of peak stresses is given in figure A3. For the large amplitude, the peak stress drops after the first cycle from 83 to 70 MPa, then the maximum value rises to 82 MPa after 6 cycles. Afterwards, there is an increase up to the 16th cycle to reach 76 MPa. For the smaller amplitude, the peak stress oscillates around 50-55 MPa.
A4. Prismatic loops: interstitial and vacancy classification Figure A4 shows the evolution of first order prismatic loops separated into vacancy and interstitial types. For the large amplitude, the curves corresponding to both types show an  In the case of the small amplitude, up to 110 cycles both types are almost equal (≈12 prismatic loops each), afterwards no increase in interstitial type prismatic loops is seen, whereas the vacancy type continues increasing to reach a maximum of 25 prismatic loops.
A5. Influence of the free surface on slip system activity During the first cycle, the dislocation-dislocation interaction is negligible and the total stress field in the tensile axis direction y is equal to the resulting stress field due to the boundary condition ŝ yy (figure A5). After one cycle, dislocations with Burgers vector b 4 are hindered by a backstress due to the stress field s yy from dislocations of Burgers vector b 5 ( figure 5(b)).

A6. Single slip setup
The experimental work of [63] was simulated for an initial dislocation density of ρ t =2×10 12 m 2 (initial source-length=4745a) using only the glide system with the  highest Schmid factor (glide system 9) as the initial dislocation microstructure and a plastic strain amplitude of ε pl =0.08%. The crystal system was chosen such that the [982]-direction was aligned along the y-axis and the [22̅ 1̅ ]-direction along the x-axis. Figure A6 shows the dislocation microstructure after 3 cycles. The simulation boundary conditions and grain geometry are the same as in figure 1. ORCID iDs D Weygand https://orcid.org/0000-0002-8681-3904