Data-driven exploration and continuum modeling of dislocation networks

The microstructural origin of strain hardening during plastic deformation in stage II deformation of face-centered cubic (fcc) metals can be attributed to the increase in dislocation density resulting in a formation of dislocation networks. Although this is a well known relation, the complexity of dislocation multiplication processes and details about the formation of dislocation networks have recently been revealed by discrete dislocation dynamics (DDD) simulations. It has been observed that dislocations, after being generated by multiplication mechanisms, show a limited expansion within their slip plane before they get trapped in the network by dislocation reactions. This mechanism involves multiple slip systems and results in a heterogeneous dislocation network, which is not reflected in most dislocation-based continuum models. We approach the continuum modeling of dislocation networks by using data science methods to provide a link between discrete dislocations and the continuum level. For this purpose, we identify relevant correlations that feed into a model for dislocation networks in a dislocation-based continuum theory of plasticity. As a key feature, the model combines the dislocation multiplication with the limitation of the travel distance of dislocations by formation of stable dislocation junctions. The effective mobility of the network is determined by a range of dislocation spacings which reproduces the scattering travel distances of generated dislocation as observed in DDD. The model is applied to a high-symmetry fcc loading case and compared to DDD simulations. The results show a physically meaningful microstructural evolution, where the generation of new dislocations by multiplication mechanisms is counteracted by a formation of a stable dislocation network. In conjunction with DDD, we observe a steady state interplay of the different mechanisms.

The microstructural origin of strain hardening during plastic deformation in stage II deformation of face-centered cubic (fcc) metals can be attributed to the increase in dislocation density resulting in a formation of dislocation networks. Although this is a well known relation, the complexity of dislocation multiplication processes and details about the formation of dislocation networks have recently been revealed by discrete dislocation dynamics (DDD) simulations. It has been observed that dislocations, after being generated by multiplication mechanisms, show a limited expansion within their slip plane before they get trapped in the network by dislocation reactions. This mechanism involves multiple slip systems and results in a heterogeneous dislocation network, which is not reflected in most dislocation-based continuum models. We approach the continuum modeling of dislocation networks by using data science methods to provide a link between discrete dislocations and the continuum level. For this purpose, we identify relevant correlations that feed into a model for dislocation networks in a dislocation-based continuum theory of plasticity. As a key feature, the model combines the dislocation multiplication with the limitation of the travel distance of dislocations by formation of stable dislocation junctions. The effective mobility of the network is determined by a range of dislocation spacings which reproduces the scattering travel distances of generated dislocation as observed in DDD. The model is applied to a high-symmetry fcc loading case and compared to DDD simulations. The results show a physically meaningful microstructural evolution, where the generation of new dislocations by Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. multiplication mechanisms is counteracted by a formation of a stable dislocation network. In conjunction with DDD, we observe a steady state interplay of the different mechanisms.
Keywords: crystal plasticity, continuum dislocation dynamics, dislocation networks, data science, data-driven modeling (Some figures may appear in colour only in the online journal)

Introduction
Dislocation networks interconnecting different slip systems are a key feature of stage II hardening in fcc single crystals which have been observed and described in early experimental works [1][2][3]. Recent discrete dislocation dynamics (DDD) investigations allow for a deeper understanding of the microstructural features which govern the effective mobility of dislocations in networks [4]. It has been shown that dislocation multiplication mechanisms lead to cascades of further reactions and multiplication events, which can span the whole simulation volume even though the multiplication events are local. In the course of this observation, a dislocation may trigger several other reactions but its own travel distance is limited by the reaction it triggers. The consequence is the formation of a dense and relatively stable dislocation network. Within the network, a large scatter in dislocation spacings and travel distances of individual dislocations is observed [4,5]. Thus, dislocation networks are characterized by individual and very heterogeneous dislocation reaction and multiplication events. Some existing models, e.g. [6,7], approach the challenge of modeling dislocation-based plasticity in the stage II hardening regime by phenomenological formulations based on modifications of the Kocks-Mecking formulation [8] using a 'Taylor'-like mobility law [9] to capture strain hardening. Other models account for immobile dislocation density, which either incorporates forest interaction as presented e.g. in [10][11][12], or the formation of dislocation dipoles, e.g. in [13,14]. However, most existing approaches combine all mechanisms into slipsystem-wise formulations, usually incorporating an internal length scale which is referred to as 'mean free path' of dislocations [10,11]. Therefore, microscopic details such as the dislocation density increase on inactive slip systems and the explicit coupling of several slip systems by dislocation reactions are not reflected. Some steps in direction of including more microstructural features have been achieved by temporal statistics of dislocation ensembles, e.g. by [15] which is applied to a continuum formulation in [16], and by including aspects of dislocation multiplication, e.g [7,17]. However, a consistent continuum formulation of dislocation network formation is still missing. Moreover, despite the usual interpretation of the 'Taylor' relation as the critical shear stress to move 'pinned' dislocation segments [18,19], statistical variation of segment lengths in dislocation networks have not received much attention in existing continuum models.
In order to reflect the dislocation dynamics in dislocation networks, we aim for a clear separation of the dislocation line length increase, resulting in plastic slip, from the mechanisms which govern the rate of the plastic slip. In contrast to the edge-screw approaches based on [20], the continuum dislocation dynamics (CDD) theory [21][22][23] formulates the kinematic properties of curved dislocation lines. This is achieved by incorporating the so-called curvature density additionally to the dislocation density, where the former represents the change in dislocation curvature along the dislocation line. Since the integral of the curvature density can be interpreted as the number of dislocations, the CDD theory allows for the explicit consideration of mechanisms which separately address the dislocation line and the dislocation as a closed object. This separation provides the basis for the formulation of dislocation multiplication by glissile reactions and cross-slip in [24], which is a first step towards a physically based formulation of dislocation networks in a dislocation based continuum theory. However, in order to reflect further details of dislocation network evolution, dislocation-based continuum models still lack relevant coupling mechanisms between slip systems which lead to stable networks. With regard to a physically based modeling, the identification of meaningful relationships in dislocation networks is needed and a concept of how to approach the statistical variance in dislocation spacings and travel distances has to be developed.
We present a model of dislocation network formation based on physically inspired formulations for the evolution of dislocation reactions while maintaining the kinematic description of dislocation expansion and plastic strain evolution based on single slip systems. Using data-driven methods for the analysis of DDD networks and systematic feature selection to identify physically meaningful correlations, further insights into the complexity of dislocation networks are obtained and continuum terms are derived to represent dislocation reactions during the network evolution. Based on the found relationships between microstructural features, we extend the model of dislocation multiplication presented in [24], by formulating evolution equations for dislocation reactions leading to sessile dislocation junctions (Lomer reaction) and junctions of zero Burgers vector (collinear reaction). The proposed continuum model accounts for two competing mechanisms: (i) the generation of new dislocations, which allows for further plasticity and (ii) the limitation of the dislocation mobility, i.e. the limitation of the plasticity generated by each individual dislocation. The interplay of the different mechanisms is analyzed in a simple system followed by a comparison of the model to DDD simulations [4]. As a result, dislocation network formation is achieved by explicit consideration of the evolution of dislocation reactions-through their formation and dissolution-rather than generalizing the mean dislocation distance as the dominant dislocation-mobility criterion for describing the properties of dislocation networks.

Voxel-based analysis of dislocation networks
The strain hardening in stage II deformation of face-centered cubic (fcc) crystals is largely determined by effects originating from the intersection of dislocation lines on different slip systems, which leads to the formation of dislocation networks. According to DDD observations in [4,5], the most important mechanisms which govern the dislocation microstructure in such networks include cross-slip as well as glissile, Lomer and collinear reactions. Based on the 3d-DDD dislocation structure investigated in [4], we consider a tensile test using a fcc single crystal with a volume of (5 μm) 3 . The system is subjected to a constant strain rate ofε = 5000 s −1 , the initial dislocation density is ρ ≈ 1.15 × 10 13 m −2 . For a detailed description of the DDD material data see [4]. With respect to a subsequent continuum consideration, we transform the DDD data, consisting of discrete dislocation lines, into voxels, which contain the continuum field variables as used in the CDD formulation according to [23].

Preparation of the voxel data set
We subdivide the simulation domain into evenly-sized cubic voxels, as visualized by figure 1. Along each spatial dimension, we use five voxels, thereby getting a voxel size of (1μm) 3 . We represent each voxel by characteristics extracted from the DDD data. These characteristics are continuum field variables, such as the total dislocation density ρ tot , or the reaction densities. The term reaction density denotes the dislocation density representing the line length of the reaction in the specific voxel, i.e., the line length of glissile (ρ gliss ), Lomer (ρ Lomer ) or collinear reactions (ρ coll ) summed over the slip systems. In the DDD framework, dislocation reactions are realized by a superposition of two segments. The resulting junction remains in the data structure, which allows to track the evolution of all junction types. We ignore Hirth reactions, since they are only very rarely observed in the underlying DDD simulations and show hardly any contribution to strain hardening, as also observed in [5]. Further, we restrict the present investigation to mechanisms in dislocation networks which are caused by intersection of dislocation lines leading to dislocation reactions. Therefore, we do not explicitly analyze dislocation mechanisms which typically dominate in single slip conditions, such as dislocation dipole formation. We also do not choose cross-slip as a voxel characteristic in the present section, since its formation mechanism is largely different from dislocation reactions. Moreover, we do not investigate the dislocation velocity in the DDD data, since its discrete characteristic in DDD is not comparable to the homogenization of the velocity in CDD in which the velocity is a spatial variable originating from averaged stress terms (see e.g. [25]). Regarding the temporal dimension, we choose sufficiently large time intervals to average local scattering in the DDD data, resulting in a number of 22 time intervals overall. For each time interval, we average over 50 DDD time steps, resulting in a difference in time of ≈ 5 × 10 −9 s between time intervals. Considering the spatial and temporal dimensions, we get 2750 data objects. However, the actual number of voxels used in our investigation is lower and depends on the reaction type, as we exclude voxels with a reaction density of zero, i.e., voxels without dislocation reactions. Depending on the reaction type, the actual number of data objects thus reduces to between 91% and 87% of the full data set.

Exploration of the data set
In this section, we investigate the evolution of the reaction density as obtained from the data set. Figure 2 shows the evolution of the measured line length of the DDD simulation stored in dislocation reactions for a discretization of five voxels along each spatial direction. We represent each voxel and time step individually, but sum over the slip systems. The plot shows a comparable increase of reaction density with the dislocation density for all reaction types. The high relative initial scatter decreases with an increasing total density. The overall increase of reaction density with dislocation density motivates our prediction models in the next section. Additionally, the coloring in figure 2 shows the number of dislocations per voxel and time interval. The value is determined by counting the number of closed dislocations with unique loop index within each voxel, such that one dislocation can be part of several voxels, but each loop index is only counted once in each voxel. It can be seen that the number of dislocations in each voxel also increases with increasing total dislocation density and reaction density.

Prediction of reaction density
Based on the findings obtained in the previous section, we want to find an adequate prediction model for the reaction densities. Considering a reaction density ρ react , representing glissile, Lomer or collinear reaction density, we assume the following functional form: where π react denotes the set of all pairs of reacting slip systems for the particular reaction type. ξ and ζ are indices of the slip systems. The function relates the line length of all reactions present at a given dislocation state to the overall dislocation content and dislocation spacing on either of the coupled slip systems. Therefore, as a natural assumption for incorporating intersecting dislocation lines on two slip planes, the dislocation densities on both slip systems are coupled with the dislocation spacing on the respective forest system, which leads to the √ ρterms. This relationship can be motivated independent of the reaction type as coupling terms of the slip system pairings π react . We assume that this functional form is applicable regardless of the type of reaction, if no change in loading path is involved. As an objective of this section, we evaluate the validity of this assumption based on data-driven analyses. The reaction pairs depend on the reaction type: we have 24 reaction pairs for ρ gliss , 12 reaction pairs for ρ Lomer and 6 reactions pairs for ρ coll , which directly originate from the fcc crystal structure. The coefficients β ξζ 1 and β ξζ 2 are unknown yet. These coefficients specify the relationship between the dislocation density on the reacting slip systems and the line length of resulting dislocation reactions.
As the equation is linear in the interaction terms, ρ ξ ρ ζ and ρ ζ ρ ξ , we validate the equation by training and evaluating a linear regression model [26] separately for each reaction type. By using multi-linear regression, we estimate the coefficients β ξζ 1 and β ξζ 2 of equation (1) with a least-squares fit. The reaction density is the prediction target and the corresponding interaction terms are the features put into the model, respectively. Before training the regression models, which we obtain from the scikit-learn library for Python, we apply min-max scaling to get all features to a common range. Furthermore, we remove features which have an absolute Pearson correlation of at least 0.95 to another feature, as linear models are known to have problems with features that are linearly related [26]. After training, we evaluate all models with the proportion of variance predicted, R 2 . To determine how well our models generalize, we apply a train-test split along the temporal axis: data from the first 17 time intervals forms the training set and data from the last 5 time intervals forms the test set. Figure 3 compares the values predicted by the trained regression models with the ground truth values. For a discretization of five voxels along each spatial direction, we observe a high prediction quality for all three reaction types, underpinning the relationship between dislocation density and reaction density. This also shows the validity of the earlier assumption to apply the same functional form independent of the specific reaction type.
To further investigate the relationship between dislocation density and reaction density, we also explore four alternatives to equation (1). Since the values of β ξζ 1 and β ξζ 2 differ even within the same reaction pair, when using equation (1), we first use identical coefficients for each two terms in a reaction pair, i.e., ∀ξ ∀ζ β ξζ 1 = β ξζ 2 , to analyze the impact of that simplification. Second, we simplify equation (1) even further by using the same coefficient for all interaction pairs: Third, we combine equation (1) with a simple feature selection technique to get a model with only a few coefficients: for all reaction types, we sort the interaction terms according to the correlation to the prediction target, i.e. reaction density, and only choose one third or one sixth of the interaction terms with the highest correlation as features. We then train the regression model with only these features, i.e. with a limited number of interaction terms. Fourth, we use the raw dislocation densities from the 12 slip systems instead of the interaction terms. Table 1 shows the results for the predictions with the different feature sets. We observe that setting β ξζ 1 = β ξζ 2 does not reduce the prediction quality. In contrast, merging all interaction terms according to equation (2) results in a drop of prediction quality. However, as a trade-off, the number of coefficients is reduced significantly, making the model simpler. The same holds for the feature selection approach, which results in a prediction accuracy between equations (1) and (2). Furthermore, we can also see an accuracy decrease when increasing the feature selectivity, i.e., selecting 1/6 of the interaction terms instead of 1/3. As a simple sanity check, we also tested selecting the same number of features, but taking the features with the lowest absolute correlation to the target instead of the highest. As one could expect, prediction accuracy turned out to be clearly worse. Finally, using the raw dislocation densities without the pre-defined interactions results in a similar or worse prediction performance compared to the feature sets using equation (1), with the largest difference for the collinear reaction density. Figure 4 displays the correlation of features from different feature sets with the target variable, i.e. each reaction density. Here, we observe that the raw dislocation densities overall show a lower correlation with their prediction target than the interaction terms from equation (1). However, one needs to be aware that these correlations only analyze the relationship between each feature on its own and the target variable. In contrast, the multiple linear regression model uses a combination of all features from the corresponding feature set. This methodological difference can explain the discrepancy between correlation results and prediction results regarding the combined interaction term of equation (2). In this case, there is only one feature per target variable and therefore only one correlation value. We observe a high correlation with the target variable, but still obtain no better prediction quality than the individual interaction terms from equation (1) combined in a model.

Dislocation networks in a dislocation-based crystal plasticity framework
3.1.1. Elasto-plastic framework. In the continuum, we describe the elasto-plastic deformation of fcc metals by an additive decomposition of the small-strain distortion tensor into an elastic and a plastic part using the gradient operator D. The plastic distortion, described by the tensor β pl , is assumed to originate solely from the evolution of the dislocation state and is determined by evaluating the sum of the plastic slip γ ξ (see section 3.1.2) over all slip systems ξ Herein, the orientations of the slip systems are defined by the orthonormal basis {d ξ , l ξ , m ξ } with the slip plane normal m ξ and the slip direction Following [21,22], the evolution of the dislocation microstructure can be described by the evolution of the total dislocation density ρ ξ , the vector of the geometrically necessary dislocation density κ ξ = (κ ξ screwd , κ ξ edge , 0), and the curvature density q ξ . The latter variable describes the local change in angular orientation of an ensemble of dislocation lines in an averaging volume. The dislocation flux-based kinematic evolution of the dislocation state-neglecting any change of dislocation content by dislocation reactions or cross-slip-is therefore described by the following equations, assuming an isotropic dislocation velocity v ξ : Here, we employ the closure assumptions introduced in [21] using the dislocation alignment tensor The relationship between shear stresses on the slip system and the dislocation velocity v ξ is characterized by a velocity law, assuming a linear dependency between velocity and effective resolved shear stress τ ξ . The effective resolved shear stress originates from a superposition of an external loading or boundary conditions projected on the slip plane and internal stresses of the dislocation microstructure: τ ξ = τ ξ ext + τ ξ int . The internal stresses represent the mean dislocation stress fields using a mean field approach as in [25] according to the numerical resolution of the system, and short range correction stresses. However, the long-range internal stress fields vanish in configurations which only consist of statistically stored dislocations (SSDs). This problem can be solved by introducing formulations for dislocation interaction and reaction, as in [14,24] and as explained later in this section.
The elasto-plastic formulation described in this section is incorporated in the finite element code M++ based on a parallel multigrid method [27,28], which further incorporates the numerical framework for the solution of the dislocation microstructure problem described by equation (5). The CDD framework is described in detail in [23]. Using a discontinuous Galerkin scheme, the equations are solved separately on each slip system. The elasto-plastic problem is solved by a standard finite element solver, where the same mesh is used for the stress computation and the evaluation of the internal variables of the microstructure.

Dislocation networks in the continuum.
Due to the planar nature of the kinematic evolution of the dislocation state in CDD (equation (5)), any mechanism which originates from the reaction of dislocations on different slip systems has to be included by rate equations, which modify the dislocation content as well as the curvature density on the respective slip systems. This is described in [24] for the aspect of dislocation multiplication. As for the voxelbased analysis of the DDD dislocation network in section 2, we formulate physically based evolution equations for dislocation reactions, which will be derived in detail in section 3.2 in the context of the given continuum theory. A schematic visualization of all considered mechanisms is shown in figure 5. Cross-slip and glissile reactions lead to dislocation multiplication by generation of new dislocations [4,24]. In contrast, Lomer reactions lead to sessile dislocation junctions which typically impede the motion of the involved dislocations [29,30]. The latter process is further called stabilization of dislocations. Collinear reactions lead to an annihilation of collinear dislocation lines sections. All reactions involve two dislocations on different slip systems, which intersect and form a dislocation reaction, as long as at least one segment is mobile.
In DDD simulations [4], it is observed that dislocations mostly do not further contribute to plastic slip once incorporated into the network e.g. by formation of dislocation reactions. Further, the network microstructure is characterized by a large scatter of the velocities and travel distances of involved dislocations. The resulting length of dislocation lines connecting endpoints of dislocation junctions has been referred to as link lengths, where statistical variations of this link length have been proposed as an explanation for the emerging microstructure of dislocation networks [5,31]. Therefore, the onset of plastic flow can be related to a 'weakestlink' argument [32]. In the context of a dislocation-based continuum theory, we therefore incorporate a formulation in which the total dislocation density ρ ξ is separated into a mobile dislocation density ρ ξ M , which evolves according to equation (5), and into a part which we call network dislocation density ρ ξ net , which does not contribute to plastic slip: For the rate of the plastic shear γ ξ it thus holds The network dislocation density ρ ξ net is distributed over a certain range of dislocation link lengths to capture the statistical variations in travel distances of the dislocations. This aspect will be derived in detail in section 3.2.2. The density ρ ξ net in equation (7) is further subdivided into (i) a density of Lomer junctions ρ ξ Lomer , which captures the line length of the junction between both end-nodes (visualized by the blue line in figure 5(c)) and (ii) a stabilized dislocation density ρ ξ S . The latter contains the line length of dislocations which are part of the stable network due to their attachment to Lomer junctions (visualized by the red and green lines in figure 5(c)). Both densities can be increased by the formation of Lomer junctions. This results in Here, the density of Lomer junctions ρ ξ Lomer is shared with the second involved slip system, resulting in a prefactor of 0.5. The unzipping of Lomer junctions is reproduced by a transformation of the respective fraction of the network dislocation density (ρ ξ net ) back into mobile dislocation density, as will be derived in section 3.2.
For the evolution of the dislocation network, the network dislocation density ρ ξ net is defined as an internal variable, which evolves based on sink and source terms in the evolution equation of the mobile dislocation density ρ ξ M . The evolution of the averaged dislocation state using equation (5) and including cross-slip, as well as glissile, Lomer and collinear reactions results in Here, () react is an index denoting different dislocation reactions, i.e. () react = () gliss , () Lomer , () coll . Further, any density increase, resulting from the formation of a reaction is denoted by( ), whereas the respective increase resulting from unzipping of Lomer junctions is denoted bŷ ( ). Dislocation reactions and cross-slip lead to a decrease of dislocation density on the reacting slip systems, which is denoted by ∂ t ρ ξ M,react , ∂ t ρ ξ M,cross , ∂ t ρ ξ S,react . Dislocation multiplication by cross-slip and glissile reactions leads to a generation of new dislocations which increases the mobile density ∂ tρ and the screw-GND density ∂ tκ ξ cross in case of cross-slip (see also [24]). The density of Lomer junctions is reflected by 0.5∂ t ρ ξ Lomer , where the reaction further transfers dislocation line length from mobile dislocation density (ρ ξ M ) to its network counterpart (ρ ξ net ) by ∂ tρ ξ M,Lomer due to the stabilization of attached dislocation lines. The reverse process is the unzipping of Lomer junctions, which increases the mobile density by ∂ tρ ξ M and equally decreases the network density by ∂ tρ ξ net = ∂ tρ ξ Lomer + ∂ tρ ξ S . All dislocation reactions and cross-slip lead to a decrease of dislocation density on the reacting slip systems ∂ t ρ ξ M,react , ∂ t ρ ξ M,cross , ∂ t ρ ξ S,react due to the mechanism itself. As derived in [24], we consider a concentration of dislocation curva- ture into the end-nodes of dislocation junctions which is reflected by ∂ t q ξ react and ∂ t q ξ cross , see also section appendix C. Table 2 provides an overview of mechanisms, which can form mobile or network dislocation density. These mechanisms are derived in detail in section 3.2.
The explicit incorporation of Lomer reactions in the dislocation network (equations (7) and (9) governs the evolution of mobile dislocation density (ρ ξ M ) and thus the plastic slip evolution defined by the Orowan equation equation (8). Restricting the increase of mobile density by junction formation therefore acts as a source of hardening. Any increase of the density of the dislocation network (equations (7) and (9)) inhibits the plastic slip and thus has an effect similar to a local yield stress. This idea is conceptually comparable to formulations which introduce dislocation dipoles as a source of hardening, e.g. in [13,14]. In addition to the strengthening due to junction formation, existing Lomer junctions (described by ξ ρ ξ Lomer ) can act as obstacles to moving dislocations due to their intrinsic stress-field. This is considered by an additional hardening component ∝ ζā Lomer ρ ζ Lomer with a Lomer interaction coefficient a Lomer contained inā Lomer = 0.5a Lomer . Here, the factor 0.5 originates from the assumption that the Lomer density is shared between two slip systems. For the Lomer coefficient we choose a Lomer = 0.122 for simplicity, which is the same value as the Lomer and self hardening coefficients in [6,19], but represents a stress interaction with mobile dislocations instead of the strength due to the junction formation. For glissile and collinear reactions the additional hardening does not exist, since in case of the collinear reaction the lines annihilate and for glissile reactions the dislocation is mobile and included in the mobile dislocation density. The contribution of glissile, collinear and Hirth reactions as well as self and coplanar interaction to hardening is considered by an interaction matrix a ξζ according to [33]. The values of the hardening components are taken from [6] where the Lomer component is set to zero and replaced by the explicit incorporation of Lomer reactions (equations (7) and (9)). In combination with the interaction of mobile dislocations with existing Lomer junctions explained above, we thus obtain with the shear modulus G. The stress term τ ξ fl, mat is incorporated into the velocity law, using the effective resolved shear stress τ ξ , as with a friction coefficient of B = 1 × 10 −4 Pa s.

Formation of dislocation reactions.
In the course of the voxel-based investigation of the DDD data in section 2, a function between dislocation density and reaction density was motivated, which couples the dislocation activity on respective slip systems, see equation (1). Here, a formulation of dislocation reactions for a continuum theory in conjunction with the investigation in section 2 is derived. Since equation (1) formulates a relation for the reaction density present at a given dislocation state, the equation can not be directly transferred into a continuum framework which relies on a formation rate of dislocation reactions. However, we follow the argumentation in section 2 regarding the validity of the investigated function for all reactions shown schematically in figure 5 to formulate reaction equations which incorporate a slip rate according to [24,34]. Regarding dislocation multiplication, a crossslip rate derived from temporal statistics of discrete dislocation ensembles has been applied to a continuum framework by [16]. In order to provide a consistent formulation of dislocation networks, we incorporate dislocation multiplication by cross-slip and glissile reaction according to the model presented in [24]. For the derivation of the dislocation reaction model in the context of the dislocation-based plasticity framework explained in section 3.1, we now consider the pairing of two reacting slip systems ξ r and ξ r individually. As a limiting criteria, at least one of the involved dislocation lines has to be mobile, i.e. ρ r M , ρ r M > 0 and v r > 0 or v r > 0. The stabilized dislocation density ρ ξ S is also a possible partner for dislocation reactions, although it does not contribute to plastic slip. Further, we assume for simplification that the statistical variation in dislocation travel distances results in a distribution of the dislocation link lengths around the mean dislocation spacing, such that the density of all link lengths, i.e. line length times number of dislocations per volume, is constant. Using those assumptions, the mean distance of the involved dislocations on the intersecting slip systems ζ is L ζ = 1/ ρ ζ M + ρ ζ S . The increase in reaction density is then described by a term comparable to [24,34]: with () react = () gliss , () Lomer , () coll . Here, C react is a combined coefficient similar to β in equation (2) which describes the effective length of the dislocation junction and contains the Burgers vector. The value range of this coefficient has been determined by DDD-simulations for glissile reactions, i.e. C react = C gliss [34,35]. Since equation (13) is valid for all reactions, the density ρ react is not tied to a specific slip system and therefore does not include a slip system index.
The reaction equation (equation (13)) describes the amount of density which is affected by the intersection of two dislocation lines forming a reaction, but it does not state if the reaction is stable. We therefore describe the stability of any junction, i.e. the irreversibility of the reaction formation, by introducing a fraction of stable dislocation reactions η ξ react for all reaction types, i.e. () react = () gliss , () Lomer , () coll . The process of a reaction formation is assumed to be irreversible for glissile and collinear reactions. Thus, it holds η ξ gliss = η ξ coll = 1. Whereas the collinear reaction leads to a dislocation annihilation, which is effectively a zero Burgers vector junction, the glissile reaction leads to the generation of a new mobile dislocation on a third involved slip system denoted by ξ gl . It thus holds ∂ t ρ gl gliss = η gl gliss ∂ t ρ gliss . Following [24], the glissile reaction leads to an increase in mobile dislocation density and in curvature density due to a bow-out of the junction up to half of the mean dislocation spacing (see appendix B for details).
In contrast to the glissile and the collinear reaction, the Lomer reaction creates a sessile dislocation junction which can unzip if the effective resolved shear stress τ ξ exceeds a critical shear stress τ cr,Lomer . The latter depends on the length of the attached dislocation segments [36,37]. We therefore assume the fraction of stable Lomer reactions η ξ Lomer to scale within the interval [0, 1]. The increase of the density of Lomer junctions is then described by using ρ Lomer from equation (13) with the index () Lomer . Although the Lomer junction crystallographically does not belong to either of the reacting slip systems, we assign the density of Lomer junctions to both involved slip systems and weight its contribution to the network dislocation density by a factor of 0.5, see equation (9).

Stability of Lomer reactions. The fraction of stable Lomer reactions η ξ
Lomer takes into account that Lomer junctions can dissolve if a large enough shear stress, i.e. a critical stress τ cr,Lomer , is applied. The strength of the junction depends on the link length, i.e. the length of the attached dislocation lines [36,37], which is assumed to scale with the dislocation density by 1/ √ ρ. It thus holds η ξ Lomer = f τ ξ , ρ ξ . Therefore, an intersection between two dislocation lines on two slip systems does not result in a stable Lomer junction if τ ξ > τ cr,Lomer . In the continuum model, η ξ Lomer then acts as a prefactor in the reaction equation equation (14). The critical stress τ cr,Lomer can be approximated to τ cr,Lomer ≈ 0.5Gb 1 L Lomer (15) based on line-tension arguments [37,38], or based on atomistic and DDD-simulations [19,36,39,40]. Here, L Lomer denotes the link length, for which we assume L Lomer ∝ 1/ √ ρ tot = 1/ ξ ρ ξ . As introduced in section 3.2.1, we assume a distribution of link lengths due to scattering travel distances of dislocations. On both reacting slip systems ξ r and ξ r , we thus find a maximum (L max ) and minimum length (L min ) of dislocation lines, which results in a distribution of the dislocation link lengths in the averaging volume between L min < L = 1/ ξ ρ ξ < L max . Figure 6(a) schematically shows a dislocation configuration for L min and L max . In the following, we choose L min = 0.1L and L max = 2L, based on the measured scattering of the travel distances and the nucleation radii of the dislocations within the DDD dislocation network [4].
Due to the relation between critical shear stress and dislocation link length (equation (15)), there exists a longest dislocation line segment which can be stable for a given resolved shear stress . τ denotes the applied shear on either slip system. The overall range is limited by the maximum (L max ) and minimum (L min ) length in conjunction with the minimum (τ min ) and maximum (τ max ) shear stress (see equation (16)). The maximum stable segment length is limited by L ξ cr,Lomer = min L r cr,Lomer , L r cr,Lomer on the coupled slip systems.
This defines the weakest Lomer reaction. Figure 6(b) schematically shows the length distribution of the dislocation segments for the general case of a different resolved shear stress on both reacting slip systems. The stress range in which stable Lomer junctions can exist is limited by L min (which defines the strongest junction) and L ξ cr,Lomer with L ξ cr,Lomer = min L r cr,Lomer , L r cr,Lomer , as visualized by the green area. This means the stability of a Lomer junction is limited by the slip system with the highest shear stress. The fraction of stable Lomer reactions is then defined by An increase in shear stress on at least one slip system therefore results in ∂ t L ξ cr,Lomer < 0. This reproduces the physical observation of unzipping Lomer junctions, which we simplify by assuming the reverse process of the Lomer reaction formation described by equation (14). This means a transformation of the fraction ∂ t L ξ cr,Lomer / L ξ cr,Lomer − L min of the density of Lomer reactions, denoted byρ ξ Lomer , and the stabilized dislocation densityρ ξ S into mobile dislocation densityρ ξ M . Thus it holds where ρ ξ Lomer and ρ ξ S is the density of Lomer reactions and the stabilized dislocation density in the given time step.

Impact of dislocation reactions on the dislocation line length.
In conjunction with [4,24], we assume that dislocation reactions as well as cross-slip lead to a loss of line length on the reacting slip system. This involves line length which vanishes in the junction, i.e. the line between both junction end-nodes. The subdivision of the total dislocation density into a mobile and a network part in equation (7) implies that the respective loss of line length has to be broken down according to which dislocations are able to react. This can be explained by an example: the mobile dislocation density on the forest slip system ρ ζ M can only reduce the dislocation densities on the primary slip system ρ ξ M and ρ ξ S by their relative amount ρ ξ M / ρ ξ M + ρ ξ S and ρ ξ S / ρ ξ M + ρ ξ S . The amount of mobile and stabilized dislocation density (denoted by the index () M,react and () S,react ) which is lost on reacting slip systems due to the formation of any dislocation reaction with a specific length is derived as This equation holds for all considered dislocation reactions, i.e. () react = () gliss , () Lomer , () coll . The sum of both contributions in equation (19) thus equals the negative production rate of the density of dislocation reactions ρ react in equation (13) multiplied with the fraction of stable dislocation reactions η ξ react .
The stabilized dislocation density ρ ξ S describes the line length of dislocations which are attached to the end-nodes of Lomer junctions, as shown in figure 6(a). The increase of ρ ξ S thus results from a stabilization of mobile dislocations involved in Lomer dislocation reactions. This contribution is denoted by ∂ tρ ξ M,Lomer and is assumed to be reversely proportional to the reduction of mobile dislocation density due to Lomer junction formation (equation (19) with index () Lomer ), leading to

Considered systems
3.3.1. Simplified system. First, we only consider glissile and Lomer reactions by using a configuration of three slip systems, as shown schematically in figure 7. We exclude all spatial gradients in the evolution equation equation (10), which results in a system of ordinary differential equations. The spatial orientation of the slip systems is not related to any real crystallographic orientation. Thereby, the system is artificial, but comes with the benefit to reduce the complexity in order to focus on the interplay of dislocation multiplication and dislocation network formation. The system is therefore fully homogeneous and serves as a benchmark configuration. Two slip systems ξ = 1 and ξ = 2 are subjected to a constant shear stress of τ 1,ext = τ 2,ext = 15 MPa. On the third slip system ξ = 3, we set τ 3,ext = 0, in order to mimic active and inactive slip systems. The initial dislocation density is chosen as ρ ξ = ρ ξ M = 1 × 10 12 m −2 with q ξ = 0 initially on each slip system. Due to the neglected transport terms it further holds κ ξ = 0. In order to determine the influence of Lomer junctions as obstacles for the movement of mobile dislocations, we investigate the impact of the additional interaction stress given by equation (11) using a Lomer = 0.122. For the simple system, the interaction matrix a ξζ is chosen to equal zero to solely focus on the influence of the Lomer reaction as obstacles. The reaction coefficient C react = 0.064 in equation (13) is chosen according to [24] for glissile and Lomer reactions. The elastic material parameters are given by the elastic modulus E = 71.3 GPa, a Poisson's ratio of ν = 0.34 and a Burgers vector of b = 0.256 nm.

Full fcc system with tensile loading.
In a second step, we use all 12 fcc systems and investigate the full model consisting of glissile, Lomer, collinear reactions and cross-slip. We aim to compare the average dislocation density evolution as well as the spatial density distribution to the DDD-simulations of [4]. We therefore apply the model to a cubic volume with edge length of (5 × 5 × 5)μm, as in [4,24]. The system is subjected to a tensile displacement on two opposite normal surfaces in the direction of the [100] crystal axis using a constant strain rate ofε = 5000 s −1 . This leads to a high symmetry multislip orientation in which eight slip systems have the same nonzero Schmid-factor (further called active) and four slip systems have a zero Schmid-factor (further called inactive). We set v = 0 on the inactive slip systems to avoid dislocation activity triggered by boundary conditions. The initial dislocation density distribution is chosen to reproduce the spatial features of the initial dislocation microstructure of the DDD simulations [4]. The latter originates from a relaxation of a random distribution of dislocation loops, shown in figure 8(a). The ensembleaverage of seven relaxed DDD microstructures is shown in figure 8(b). Here, the dislocation line length is averaged over 20 × 20 × 1 voxel, averaging over the full system in tensile direction (here normal to the drawing plane). It can be seen that the dislocation density accumulates in the center of the system, which results from the use of free surfaces used in the DDD simulations. For the continuum simulation, we mimic this initial microstructure by distributing the total dislocation density ρ ξ on each slip system ξ according to a multivariate normal distribution using a vector of the mean value equal to 0μm and standard deviations of σ x = σ y = σ z = 1.2 × 10 12 m −2 in the respective coordinate directions x, y and z. This results in an accumulation of dislocation density in the center of the system, as shown in figure 8(c). This initial dislocation configuration consists of network dislocation density ρ ξ net , which does not generate a long-range stress-field and therefore reproduces the relaxed dislocation network observed in [4]. For the average initial dislocation density per slip system ξ it holds ρ ξ Lomer ≈ 0.35 × 10 12 m −2 and ρ ξ S ≈ 0.75 × 10 12 m −2 , i.e. one third of the line length per slip system is stored in Lomer junctions, which resembles the configuration in [4]. This results in a total dislocation density of ρ tot = ξ ρ ξ ≈ 1.13 × 10 13 m −2 , when averaged over all slip systems. Further, we choose ρ ξ M = 0, κ ξ = 0 and q ξ = 0 initially. For the analysis, we use the interaction stress according to equation (11) with a non-zero interaction matrix in order to reproduce the mutual obstruction of coupled dislocation segments.
Following the comparison with the DDD simulations, we investigate the interplay between dislocation fluxes and dislocation network evolution by first excluding all spatial gradients (local simulation) and afterwards allowing for dislocation fluxes (nonlocal simulation). In the latter case, mobile dislocation density can leave the simulation volume through all surfaces similar to the DDD simulations. The reaction coefficients in equation (13) are chosen as C gliss = C Lomer = C coll = 0.032 for the local simulation, i.e. all reactions are weighted equally. For the nonlocal simulation, we choose C gliss = C Lomer = 0.064, C coll = 0.032 and C q = 1 in equation (C.1), which means that we increase the effect of the dislocation multiplication in order to compensate for the dislocation flux. The dislocation flux necessarily leads to a loss of dislocation density through the outer system boundaries and thus a decreasing density by using the same parameters as for the local simulation. The elastic constants are chosen as before.

Results
In order to analyze the behavior of the model derived in section 3 regarding the interplay between dislocation multiplication and stabilizing mechanisms, we first consider the simplified system explained in section 3.3.1. Afterwards, we apply the full model to the fcc system setup explained in section 3.3.2 and compare the results to DDD simulations and investigate the influence of the flux terms in the evolution equation of the dislocation system on the spacial dislocation density distribution.

Simplified system of glissile and Lomer reactions
In a first analysis, we reduce the model to only account for glissile and Lomer reactions in a configuration which consists of three slip systems, as shown in figure 7. The system is fully homogeneous since spatial derivatives in the evolution equations of the dislocation densities (equation (10)) are excluded. In a first analysis, we assume a constant dislocation  (7) and (9)) per slip system without an additional flow-stress term showing the full simulation time (a) and the initial phase of the simulation (b). (c) same as (a), but with an additional flow-stress term, see equation (11). velocity which originates from an external shear stress of τ ξ ext = 15 MPa on the slip systems ξ = 1 and ξ = 2, i.e. it holds v ξ = b B τ ξ ext . Slip system 3 is inactive, making these dislocations passive reaction partners.
The evolution of the dislocation densities representing the line lengths within dislocation networks (total dislocation density ρ ξ tot , mobile dislocation density ρ ξ M and network dislocation density ρ ξ net ) is shown in figure 9(a) per slip system over the simulation time. The total dislocation density on both active slip systems increases tenfold before showing a saturation. On the inactive slip system the dislocation density increases due to glissile reactions and reaches about half the density on active slip systems at the end of the simulation. At this point, the total dislocation density almost completely consists of network dislocation density ρ ξ net , i.e. all dislocations contribute to configurations stabilized by Lomer junctions. A close-up on the initial part of the simulation in figure 9(b) shows an initial increase in mobile dislocation density on both active and inactive slip systems. After a simulation time of about 0.4μs, the mobile density decreases along with an increase of the network dislocation density ρ ξ net . Since in this simplified configuration the system behavior is largely determined by the Lomer reaction, it is feasible to account for an obstacle-effect of Lomer junctions on remaining mobile dislocations. As explained in section 3.1.2, this is incorporated into the model by a short-range interaction stress according to equation (11) using a coefficient of a Lomer = 0.122. By choosing a ξζ = 0 in equation (11), we solely focus on the influence of Lomer junctions for simplicity. In this case, the dislocation velocity is then a function of the dislocation density, as described by equation (12). The respective evolution of the dislocation densities over simulation time is shown in figure 9(c). The system shows a relaxation into a stable network configuration, which does not qualitatively differ from the system behavior without an additional interaction stress ( figure 9(a)). However, the mobile dislocation density decreases only gradually until the end of the simulation.
The evolution of the curvature density on active slip systems displayed in figure 10(a) shows an initial increase due to glissile reactions followed by a decrease nearly until disappearance. Similar to the evolution of the mobile dislocation density in figure 9(c), the decrease of the curvature density is delayed by using the additional interaction stress. Figure 10(b) shows the evolution of the accumulated plastic slip for the simulations with and without the interaction stress equation (11). Here, a saturation of the plastic slip on active slip systems is observed after an initial increase. By using the additional interaction stress, a lower accumulated plastic slip is observed almost for the entire simulation time. However, the initial increase of the plastic slip on active slip systems is almost the same for both cases, since the initial configuration consists of mobile dislocations only.

Full fcc system with tensile loading
In a second analysis, we evaluate the performance of the continuum model including all dislocation mechanisms as visualized in figure 5 by comparing the evolution of the dislocation microstructure with DDD results published in [4]. Therefore, we apply the continuum model to the system configuration explained in section 3.3.2. For a first investigation, we exclude spatial derivatives in the evolution equation of the dislocation densities.
The evolution of the averaged total dislocation density per slip system compared to the DDD results is shown in figure 11(a) for active and inactive slip systems. The dislocation density on active slip systems increases nearly linearly after obtaining a stable plastic flow in both CDD and DDD. On inactive slip systems, the increase in dislocation density is much less pronounced as on active slip systems which matches the DDD results. Overall, a close accordance between CDD and DDD can be observed after about 0.1% strain.
The initial phase of the simulation is characterized by a mobilization of the dislocation network in both CDD and DDD. However, in the CDD simulation this leads to a short increase in density due to the dissolution of Lomer junctions, which is not observed in DDD the same way. In order to closer investigate the interplay of dislocation multiplication and the stabilizing processes in the continuum model, the breakdown of the total dislocation density into its contributions (ρ M , ρ net and ρ S ) summed over all slip systems is shown figure 11(b). The figure shows the evolution of the respective densities as fraction of the total dislocation density over strain. Here, it is clearly observed that initially no mobile dislocation density is present. Upon loading, the fraction of network dislocation density decreases rapidly to about 40% which increases mobile density. During this process, the fraction of Lomer junctions decreases from about 30% to just below 20% which explains the initial rise in density observed in figure 11(a). However, after this initial mobilization phase, the distribution between mobile and network density remains relatively stable with a division of around 60% to 40%. Corresponding results by considering dislocation fluxes in the evolution equations, i.e. the nonlocal simulation, using model parameters as outlined in section 3.3.2 show a comparable behavior, as shown in appendix A. Figure 11. (a) Evolution of the dislocation density per slip system using all explicitly modeled dislocation mechanisms, as shown in figure 5, compared to one DDD simulation and an average over 10 DDD simulations from [4]. (b) Total density subdivided into contributions from the mobile dislocation density and the network dislocation density as a fraction of the total dislocation density as a function of total strain, summed over all slip systems. Local simulation containing no dislocation fluxes.
Finally, we compare the evolving dislocation microstructure in the continuum simulation to the DDD prediction. Figure 12 displays snapshots taken from various stages of the DDD and CDD simulations, showing the spatial dislocation distributions. For a reasonable comparison between CDD and DDD, we transform the discrete dislocation lines of seven DDD simulations (shown in figure 12(a) for one simulation) into a two-dimensional representation of a dislocation density field as explained in section 3.3.2. The dislocation density distribution obtained by this averaging is shown in figure 12(b). It can be observed that dislocations which initially are already concentrated in the center of the system continue to accumulate during loading until 0.8% tensile strain. In contrast, the areas close to the boundaries remain relatively dislocation free. This applies especially to the corners and edges of the simulation volume. A similar accumulation of dislocation density is also observed in the corresponding representations of the CDD simulations, shown in figure 12(c) for the local simulation and in figure 12(d) for the nonlocal simulation. However, certain differences occur between the local and the nonlocal simulation in the later stages of the simulation. The dislocation density in the local simulation also increases in areas close to the boundaries, which is not observed in this form in either the DDD or the nonlocal CDD formulation.

Discussion
We introduce a dislocation-density-based continuum model for the evolution of dislocation networks, which combines dislocation multiplication with a formation of stable dislocation network structures. In order to find links between DDD and the continuum scale, we use methods of data science to characterize the discrete dislocation microstructure. This preliminary investigation serves as a basis for the continuum model of dislocation network evolution, where the latter is investigated in a simplified system to isolate the mechanisms, complemented by a tensile test of the full fcc system.
Given the complexity of the DDD network, manageable data sets for the analysis of the dislocation microstructure by methods of data science have to be provided. Here, the data Figure 12. Overview of dislocation microstructures obtained from DDD simulations [4] and from the continuum model for 0.0%, 0.2%, 0.4%, 0.6% and 0.8% total strain. Top view of one representative DDD simulation consisting of discrete lines (a); ensembleaverage of seven DDD simulations, averaged over 20 × 20 × 1 voxels (b) and distribution of the total dislocation density summed over all slip systems averaged over the tensile direction for the local (c) and the nonlocal continuum simulation (d). The voxelplots in (c) and (d) serve as a comparison to the averaged DDD distribution (b), but do not show the actual mesh used in the simulation, which is tetragonal. set is obtained by generating continuum field variables by averaging the dislocation characteristics such as the line length over the system volume, which is subdivided into voxels. The increase in the number of dislocations in conjunction with the reaction density on trend, shown in figure 2, indicates that the increase in reaction density mostly originates from the formation of new dislocation reactions, instead of an elongation of the junction lines. This observation motivates the assumption of a function describing the present density of dislocation reactions which is valid for all reaction types (equations (1) and (2)). By applying different scaling coefficients, i.e. β, β ξζ or β ξζ 1 resp. β ξζ 2 , we vary the aggregation of the slip system pairings and thus compare different degrees of homogenization of the dislocation interaction for the same reaction type. This is conceptually analogous to the investigation in [41] which compares the different homogenisation originating from the popular flow-stress term according to [8,9] in contrast to proposed interaction matrix by [33]. The results of the prediction with the multiple linear regression model on equation (1) demonstrate not only the existence of a relation between dislocation density and reaction density, but also confirms the validity of the assumption to use the same relationship for all reactions, as shown in figure 3. Further, the temporal split between test and training set implies that the model predicts the reaction density of future time steps based on a training on earlier stages of the simulation. The very minimal difference between the prediction on the training set compared to the test set, as shown in table 1, underlines the solidity of the model and thus the validity of the approach.
By comparing different feature sets in table 1, basically no influence of coalescing each slip system pair (equation (1) using β ξζ 1 = β ξζ 2 ) on the prediction quality is observed. For the continuum theory, this means that combining the original and reverse order of primary and forest system into a single reaction equation as done in equation (13) is a valid simplified approach. Further, the complexity of the model is significantly reduced by combining all slip system pairings of the specific reaction type into a single feature set (equation (2)). However, the drop in prediction quality indicates that a common factor independent of the pairing might oversimplify the individual contribution of each slip system. This argument is supported by the correlation between the individual dislocation densities and the overall reaction densities, which shows a very high scatter and a lower average correlation compared to the respective correlations of the coupling terms (equations (1) and (2)), see figure 4. These results indicate that while the combination of all slip systems might lead to a reasonable prediction on average, the individual contribution of each slip system and slip system pairing on the overall reaction density can largely differ. Physically this means that the line length of a forming dislocation junction is not the same for all slip system combinations. It is noteworthy that this scatter is already observed in the present high symmetry configuration which is characterized by a clear separation between active and inactive slip systems with respect to the Schmid-factors. Since the investigated setup is already the most ideal scenario of multislip loading, it can be expected that the differing influence of slip systems will increase in basically any other loading condition and potentially even in the given one for higher strains.
It can be concluded, that the results from the preliminary investigation of the DDD dislocation network further supports conclusions made earlier in that the contribution of individual slip systems to dislocation multiplication and dislocation network formation are an important feature to include in slip-system-based continuum models. What is even more significant, however, is that the in-depth analysis of the DDD data motivate similar functional relations for slip systems coupling to describe the contribution of all relevant dislocation reactions to network formation by continuum field variables. The interplay between slip systems is thus not limited to the aspect of dislocation multiplication by glissile reactions, as shown in [24,34], but also for the formation of stable structures and for collinear annihilation in dislocation networks. Since the examined functions in section 2 only relate the dislocation state to the present dislocation reaction densities, i.e. the line lengths of the junctions, the results of the voxel analysis can not be directly transferred into rate equations for the continuum framework. However, the results confirm the importance of incorporating the interplay of slip systems to reflect dislocation reactions. Therefore, we use this con-clusion as a motivation to formulate dislocation reaction equations in the CDD framework based on coupling terms similar to [24,34] which are valid regardless of the specific reaction type.
The investigation of the continuum model in the simplified system combines the aspect of dislocation expansion, which follows dislocation multiplication incorporated from [24], following dislocation multiplication with a formation of stable dislocation structures in a minimal benchmark configuration. Those mechanisms are reproduced in terms of internal variables by modifying the dislocation density due to dislocation reactions and the curvature density. For the simple system, the coupling of slip systems by dislocation reactions results in an increase of dislocation density by dislocation multiplication and a subsequent saturation due to formation of dislocation networks as shown in figure 9(a). The density increase implies a reduction in dislocation spacings, respectively link lengths, in the network, see equations (16) and (17). The saturation of the dislocation density over time is in conjunction with the saturation of the plastic slip (see figure 10(b)) caused by an increase in the fraction of stable Lomer junctions (equation (17)) at a given stress. Thereby, the Lomer junction formation upon generation of new dislocations through glissile reactions and a subsequent increase in density leads to a reduction of the link lengths (see equations (15)-(17)), limiting the plastic slip generated by each dislocation. This mechanism has been described in [4] in the context of dislocation networks and is adequately reproduced by the continuum model. Although the evolution of the total density resembles the characteristics of the similar system investigated in [24], the observed mechanism in the current study is solely caused by the reduction of the link lengths due to junction formation and is therefore the main source of hardening in the model.
The decrease of the curvature density along with the formation of the dislocation network ( figure 10(a)) reproduces the observed microstructure consisting of straight lines, which implies a concentration of the system curvature into the endpoints of dislocation junctions. This concentrated curvature is relatively stable and does not lead to an expansion of dislocations in the way it is described by the kinematic formulation of the CDD theory (equation (5)). The anti-proportional behavior of the decrease of curvature density (figure 10(a)) compared to the formation of the dislocation network (figure 9(a)) can then be interpreted as a mechanism which counteracts the increase of curvature due to generation of new dislocations and thus contributes to the stability of the dislocation network. As conclusion for the simplified system, an accordance of the proposed continuum model with observations in DDD simulations can be drawn: (I) expansion of individual dislocations is limited by their involvement in dislocation reactions. (II) The emerging dislocation network is relatively stable and consists of straight dislocation lines with varying segment lengths and concentrated curvature in dislocation junctions.
The formulation of a dislocation density which is distributed into a mobile part (ρ M ) and a part which does not directly contribute to plastic slip (ρ net ) allows for the formulation of an initial condition which reproduces stable dislocation networks. This is achieved by prescribing an initial dislocation distribution consisting solely of network dislocation density (ρ net ) in the full fcc system, shown in figure 8. Here, mobile dislocation density is not present initially and can thus only be generated by mobilization of the initial dislocation network upon loading. This explains the characteristic of the dislocation density evolution in the initial simulation phase, shown in figure 11(a). The onset of plasticity is determined by the activation of the longest dislocation lines attached to the Lomer junctions, denoted as link length, resulting in an unzipping of the junctions. Therefore, yielding originates from the mobilization of the weakest junctions, determined by the 'strength distribution' of Lomer junctions. Since the mobilization of the Lomer junction as considered here typically increases the line length, the initial yield behavior avoids a decrease of the dislocation density below the initial value in the nonlocal simulation (shown in figure A1(a) in appendix A). For simplicity, the formation and dissolution of dislocation reactions is concentrated into a single numerical time step. However, reaction processes can take place over a longer period of time and are influenced by the surrounding microstructure. This might serve as an explanation for the short increase of the density right after overcoming the yield point in CDD (figures 11(a) and A1(a)). It is also likely that a more accurate distribution of link lengths, as shown in [5], will influence the dislocation behavior in the initial part of the simulation.
Besides, an increase of the total dislocation density per slip system is observed which is very close to the DDD prediction on active and inactive slip systems for the vast majority of the simulation time for the local simulation ( figure 11(a)). Furthermore, the uniform fractional distribution of the average mobile and network dislocation densities shown in figure 11(b) indicates a stable interplay of all considered mechanisms. Thus, the dislocation multiplication is in balance with mechanisms that either annihilate (collinear reaction) or stabilize line length (Lomer reaction) within the network. The results show that this interplay of dislocation reactions is able to reproduce observations in DDD simulations [4], where it is shown that cascades of multiplication events produces further plasticity, which in turn is limited by other reactions. This conclusion holds also for the nonlocal simulation shown in figure A1(b) in appendix A. However, the average dislocation density increase is slightly nonlinear, as shown in figure A1(a). The system areas closer to boundaries contain relatively low density (see figure 12(d)), but carry a large percentage of the total volume. Therefore, even a slight overestimation of the dislocation activity in these areas can have a large impact on the total density. This overestimation is likely to originate from an inaccurate reproduction of dislocation mechanisms at very low dislocation densities in the continuum.
This difference in average density evolution motivates a deeper examination of the spatial dislocation density distribution obtained from the CDD simulations, which is shown in figure 12. The comparison with averaged DDD results demonstrates that the local and the nonlocal CDD simulation is able to resolve the concentration of the dislocation network in the center of the system and its continuous densification. However, the local simulation (figure 12(c)) predicts a noticeable density increase at the boundaries of the system, which is not observed in either the DDD or the nonlocal CDD simulation. In this respect, the distribution obtained from the nonlocal simulation (figure 12(d)) is closer to DDD due to the outflow of the density at the boundaries and a transport of dislocations into the center of the system. It might not be very surprising that the nonlocal simulation as opposed to the local simulation is able to capture such effects, since they ultimately arise from the mere existence of the nonlocal terms. However, it should be noted that the examined simulation volume with (5 μm) 3 is relatively small for a continuum theory. It is naturally to be expected that any boundary effect will have lower impact in larger simulation volumes. In conjunction with the very promising per-slip-system results for the local simulation ( figure 11) the results raise the question to what extend the incorporation of dislocation transport terms decisively determine the microstructural evolution. However, this question can not conclusively be answered in the present study.
Summarizing, the study presents an approach that combines the dislocation multiplication with the formation and evolution of dislocation networks in a continuum theory. This approach mimics the complexity of dislocation microstructures observed in DDD [4,5] and resembles similar models of dislocation dipoles in single slip [13,14]. The investigation of the averaged DDD data by data science methods shows that even in an ideal system dislocation reactions which determine the overall system behavior are influenced by individual slip systems. This opposes approaches for dislocation multiplication based on the Kocks-Mecking theory [8], as well a slip-system-wise adoption of 'flow-stress' terms in the velocity law (equation (12)) as the only source of strain hardening. Contrary to the presented model, such terms either allow for an expansion of all dislocations or stop the plastic slip completely. The main feature of the presented model is thus the separation of the mechanism that causes the plastic slip, i.e. the dislocation expansion, from the mechanisms that modify the rate of the plastic slip either by generating new dislocations or by limiting their expansion on average. Due to the explicit incorporation of dislocation reactions using their crystallographic slip system pairings, dislocations on all slip systems have a significant impact on the mobility of dislocations on individual slip systems. Furthermore, the effective mobility of the dislocation network is determined by a range of critical stresses to unzip dislocation junctions, which weakens the average dislocation spacing as a descriptive variable to some extend.
However, several open questions remain with respect to the analysis of dislocation networks as well as for modeling those in the CDD formulation. First, the presented model contains some assumptions and parameters, which originate from limitations of the purely planar CDD kinematic. This draws attention to the long-standing challenge of formulating a thermodynamically consistent multi-slip CDD. Further, it was assumed for simplification purposes that the reaction parameters are independent of the specific slip system pairing. However, according to the results of the DDD data analysis, different reaction parameters for different slip system pairings might be reasonable, indicating an analogy to the 'interaction matrix' compared to the 'Taylor' flow-stress. Finally, the boundary between mobile dislocations and sessile objects is typically less clear as assumed by the Lomer junctions in the current model. In this respect, it is to be determined how the character of dislocation networks is reflected in other system configurations, since real deformation processes rarely occur as homogeneously as assumed here.

Conclusion
We introduce a model for the evolution of dislocation networks in a dislocation-based formulation of crystal plasticity, which provides an approach for a homogenization of the interaction between dislocation multiplication and the stabilization of the emerging dislocation network. By using data-driven methods, we analyze a DDD dislocation network to identify microstructural characteristics and derive continuum equations for relevant dislocation reactions, which are incorporated into the continuum model using their crystallographic slip system pairings.
The analysis of the DDD network by data science methods shows that despite the ideal high-symmetry crystal orientation, the overall evolution of dislocation networks and thus the system behavior can be determined by individual slip systems and slip system pairings. This observation challenges existing models which base dislocation multiplication and strain hardening due to formation of dislocation networks on slip-system-wise approaches in which other slip systems only contribute to a single dislocation spacing. In contrast, the proposed model is based on the coupling of slip systems by dislocation reactions, resulting in an interplay between dislocation multiplication and limitation of dislocation expansion. This leads to an average density evolution as well as a local density distribution which is very close to DDD simulations. Here, mobile dislocations are generated due to cross-slip or glissile reactions, which leads to an increase in line length and thus decreases the dislocation spacings in the network. The consequence is a formation of a dislocation network consisting of stable dislocation structures with a range of critical stresses to unzip sessile junctions. Figure A2. Comparison of the dislocation density per slip system relative to the initial density averaged over the full system (solid lines from figure 1(a) normalized with the initial density) and over a central cubic volume with half of the edge length of the full system (dotted lines). of the full system. This inner volume cuts off the system areas adjacent to open borders and captures about 40% of the Gaussian initial dislocation distribution despite having a significantly smaller volume fraction. Figure A2 shows the average density evolution in this inner system area compared to the CDD results in figure A1(a) relative to the respective initial densities. Compared to the average over the full system, the relative density increase in the system center is almost linear. This confirms that the observed nonlinearity is caused by density in system areas outside the center, i.e. in areas which contain relatively low dislocation density, but still contribute significantly to the total density due to the high volume fraction. In the underlying DDD simulations, those areas are typically characterized by single, or very few dislocation lines which are likely to behave differently from dislocations in the dense center. However, the continuum model assumes that even at very low dislocation densities a fraction of the density is always involved in reactions, which probably does not sufficiency reproduce the statistical effects at the boundaries.

Appendix B. Glissile reactions in dislocation networks
The dislocation generated on the new slip system by the glissile reaction (see equation (13)) is can bow out under stress increasing both dislocation density and curvature density. This process is called 'dislocation multiplication', since it generates a new mobile dislocation instead of just increasing the line length of an existing dislocation. Following [24], the increase in mobile dislocation density and curvature density is described by: Here, we assume that the glissile reaction is followed by a bow-out of the junction up to half of the mean dislocation spacing, denoted by L bow , leading to the prefactor π/2 in equation (B.1).

Appendix C. Impact of dislocation reactions on the curvature density
In conjunction with the subdivision of the total dislocation density (equation (7)), we also assume that a related part of the curvature density does not contribute to the expansion of existing dislocations. This reproduces the observation made in DDD simulations that dislocation networks mainly consist of straight lines, where the dislocation curvature is concentrated in end-nodes of dislocation junctions [4]. Following the argumentation in [24], we reproduce this behavior by reducing the curvature density q ξ based on the argument of a reduction of 'potential' of line length increase by dislocation expansion due to the dislocation reaction. Therefore, the reduction of line length (equation (19)) necessarily leads to a proportional reduction in curvature density. We obtain the final formulation for the reduction of curvature density on the reacting slip systems ξ r and ξ r by a modification of the underlying equation in [24]: where we choose C q = 2 according to [24] unless defined otherwise.