Microstructural Comparison of the Kinematics of Discrete and Continuum Dislocations Models

The Continuum Dislocation Dynamics (CDD) theory and the Discrete Dislocation Dynamics (DDD) method are compared based on concise mathematical formulations of the coarse graining of discrete data. A numerical tool for converting from a discrete to a continuum representation of a given dislocation configuration is developed, which allows to directly compare both simulation approaches based on continuum quantities (e.g. scalar density, geometrically necessary densities, mean curvature). Investigating the evolution of selected dislocation configurations within analytically given velocity fields for both DDD and CDD reveals that CDD contains a surprising number of important microstructural details.


Introduction
Predicting and understanding the collective evolution of dislocation ensembles together with the resulting mechanical properties of crystalline materials is a long-standing goal of microstructure-based plasticity. From a computational perspective, reducing dislocation mechanics to a system of partial differential equations amenable to solution by standard continuum computational methods would be highly desirable. However, the nature of those curved, line-like defects makes the development of a continuum theory of dislocations an extremely difficult challenge [1], which today is still not fully mastered. One way of enriching macroscopic continuum descriptions with additional details of smaller length scales (e.g. with details of dislocations) are concurrent as well as hierarchical multi-scale approaches for bridging the atomistic to the continuum domains [2][3][4]. These methods aim at establishing a direct link between atomistic details of dislocations and macroscopic mechanical properties of materials. The level of detail in these approaches is very high, and length scales can be reached that are not accessible to one-scale methods as e.g. atomistic simulations. Multi-scale approaches, however, are still computationally very expensive and the information transfer between methods is numerically challenging. An alternative description of dislocation ensembles is utilized in mesoscopic methods such as the discrete dislocation dynamics method (DDD) [5][6][7][8][9][10][11][12][13][14][15]. In DDD the atomic scale is not resolved, and instead dislocations are represented as discrete, linear defects interacting according to the elastic theory of dislocations [16][17][18][19][20]. Because each dislocation is resolved individually, DDD is able to accurately describe the motion and interaction of dislocations in a very detailed manner. However, the computational cost of DDD scales with the number of dislocations considered and is therefore computationally expensive when it comes to large numbers of interacting dislocations (large referring e.g. to a density of ≥10 13 m −2 in a volume of ≥(10µm) 3 ). To overcome these limitations, one might seek continuous density-based descriptions as alternative to the representation of dislocations as discrete objects. Reducing dislocation microstructure to distinct "types" of e.g. geometrically necessary and statistically stored dislocations (GNDs and SSDs) allows strain gradient based continuum methods [21][22][23] to become independent of the number of interacting dislocations. These approaches result in a substantial gain in terms of tractable length and time scales which makes them interesting for engineering applications. The main drawback is that only few, strongly averaged aspects of systems of dislocations are preserved, and e.g. fluxes of dislocations cannot be accounted for at all. Recently emerging continuum models of dislocation dynamics (CDD) bridge the gap between these last two approaches by combining the efficiency of a continuum theory with the physically sound basis of a dislocation based description: 'Kröner-Nye tensor'-related models [24][25][26] are based on fluxes of GNDs, and continuum screw-edge representations are a relatively coarse but efficient way of representing dislocation loops [27][28][29], while models that consider densities of curved lines contain considerably more information [30][31][32][33][34][35] and recently even have been coupled to gradient plasticity models [36]. Among others, applications with high densities and/or high accumulated plastic strains particularly benefit from continuum dislocation dynamics models (e.g. [37][38][39][40][41]).
Because in a CDD-type continuum framework details of individual dislocations are accounted for in an average sense, one may question the physical validity and the predictive capability of such an approach. In principle it should be possible to benchmark CDD theories through comparison with averaged data extracted from DDD simulations. Up to date, contributions into this direction have been presented by A. El-Azab and co-workers [42,43], who investigated ensemble averages and statistical properties of dislocation microstructure e.g. in terms of dislocation density and orientation distributions. The objective of the present work is to compare CDD and DDD directly by analyzing dislocation microstructures that develop in time from equivalent initial conditions. Although such a detailed comparison is essential to validate assumptions used to derive the particular CDD formulation, no previous attempts into this direction exist (to the best of the authors' knowledge). One of the reasons for the lack of such direct comparisons may reside in the formulation of a CDD theory itself. In order to be comparable to DDD simulations in an average sense a CDD framework must fulfill a number of criteria: (i) its density measures must be derived by well defined statistical averaging steps from systems of discrete dislocations (i.e. in a bottom-up approach); (ii) its evolution equations should map an initial set of density measures onto a set of temporally evolved densities; (iii) it must be kinematically and (iv) dynamically consistent.
• Kinematic consistency describes the ability to evolve dislocations as curved and connected lines with the concomitant line length change during motion in a given velocity field.
• Dynamic consistency, on the other hand, describes the ability to predict the actual evolution of a dislocation system under externally applied mechanical load as well as mutual interactions between dislocations.
It is generally not recognized in the literature that kinematic consistency is a necessary prerequisite for dynamic consistency and considers only geometrical aspects of averaged systems of dislocations. While the problem of dynamic consistency is to date still far from being solved, the recently developed (so-called 'simplified' or 'integrated') CDD theory by T. Hochrainer and co-workers [34,35,44,45] was rigorously derived to solve the problem of kinematic consistency and has already demonstrated its applicability in a number of problems [34,35], including dislocation patterning [41]. The goal of this work is to establish a methodology for a direct and detailed kinematic comparison between DDD and CDD models and to demonstrate the accuracy of CDD through a number of benchmark problems.
The outline of this paper is as follows: in section 2 we give a brief overview of the general concept for validating continuum dislocation microstructure with discrete data used in this paper. In section 3 we introduce the mathematical foundations for the geometrical representation of curved dislocations lines and for obtaining averaged continuum quantities in an analytical as well as in a numerical way. We then briefly introduce the evolution equations for CDD and how CDD data can be numerically obtained from DDD in section 4. Section 5 investigates numerically the quality of the continuum model during time evolution of DDD and CDD microstructures for three different benchmark systems. Finally, in section 6 we demonstrate how the DDD model can be used for "data mining" in order to estimate the quality of one of the assumptions on which the CDD version used in this paper is based.

Outline of the CDD-DDD Validation Approach
The fundamental idea for validating the microstructure evolution within CDD is sketched in Fig. 1. Consider an ensemble of dislocation loops with the same Burgers vector b, occupying a domain of size L x × L y × z (Fig. 1A), where the slip system normal points in z direction. Each dislocation is located on a different slip plane and moves -driven by an external stress field -by glide only; no annihilation, cross slip etc are considered. Averaging over the height of a sub-volume z we obtain Fig. 1B, which serves as initial dislocation microstructure for our investigations. In the following sections a systematic method for extracting average fields (e.g. dislocation densities, Nye tensor, curvature, etc) from systems of discrete loops is developed. This discreteto-continuum conversion method (D2C ) is applied to the DDD initial configuration (Fig. 1B), from which we obtain average fields (Fig. 1C). Those serve as initial values for the CDD model. CDD initial structure (Fig. 1C) is evolved by time integration of the CDD evolution equations, while the initial configuration of discrete loops is independently evolved in time by the DDD equation of motion. This results in the DDD and CDD microstructures as shown in Fig. 1D and E. At any subsequent point of time, relevant state quantities of CDD, Fig. 1E, can be compared to corresponding DDD quantities as obtained from the D2C conversion, Fig. 1F.

Conversion of discrete dislocation lines into continuous fields (D2C)
In what follows the discrete-to-continuum (D2C ) approach is described which is used to systematically average the geometric properties of discrete dislocation lines to obtain CDD state quantities. To differentiate between variables for discrete and continuous measures the superscript 'd' is used for discrete quantities. The superscript is dropped for a continuous quantity or when the context is clear. This section starts by introducing the geometrical representation of discrete lines and their approximation through splines, followed by a formal density definition. Subsequently, mathematical averaging operators are introduced and concisely derived in a discretized form which is suitable for numerical implementation.

Discrete dislocation lines and their geometrical approximation
From a geometrical viewpoint, a dislocation line is an oriented curve c( ) which can be parameterized by its arc-length . The local unit tangent vectors are the first derivative of c while the curvature vector is given by the change of orientation per unit arc-length, The scalar curvature is the reciprocal of the local radius of curvature R of the dislocation line: Although parameterization by arc-length is the most natural choice for defining geometrical properties of a curve, in general the parameter is only available after the curve has been constructed. Thus, for practical purposes, it is often necessary to prescribe a curve through a parameter u, say u ∈ [0, 1]. In this case unit tangent and curvature become where the scalar J and its derivative are given by In this work dislocation lines are discretized into segments, each of which is represented by a cubic Hermite spline of the form where (P 0 , T 0 ) are position and parametric tangent vectors of the first knot of the segment (u = 0), (P 1 , T 1 ) are the corresponding quantities for the second knot of the segment (u = 1), γ = P 1 − P 0 α , and α is a tension parameter ‡. Cubic splines are used instead of linear segments in order to have non-vanishing curvature along each segment. This type of representation is adopted in the so-called Parametric Dislocation Dynamics method (PDD) [11]. The numerical studies below are based on the existing implementation [14,46]. ‡ Parameterization corresponding to α = 0, α = 0.5, and α = 1 are called uniform, centripetal, and chordal, respectively.

Averaging dislocation fields
As an auxiliary functional for averaging geometrical properties of discrete dislocations, we introduce the Dirac delta-operator on a dislocation line c as where the symbol '•' is a placeholder for an arbitrary line property described by a scalar, vectorial or tensorial function of c. Using (7), for each discrete line c we can define the scalar dislocation density ρ c , the line direction density (or vector of signed GND densities) c , the Kröner-Nye tensor α c [47], and the curvature density vector q c as § The later on used scalar, signed curvature density q c can be obtained by projecting the curvature vector on the outwards pointing normal where n is the slip plane normal. These discrete measures are well suited for spatial (or ensemble) averaging over a number of lines, and we introduce an averaging operator for the volume V r of size V as where the averaging volume is centered around r. Application of (13) to (8)- (12) leads to the definition of the corresponding average fields: 14) § In the following we assume only one slip system, i.e. all dislocations have the same Burgers vector b. From this case we can then construct the multi-slip system situation by appropriate superposition (e.g. summation of densities).
Therein, L V,r c ⊂ c denotes a section of line c contained inside the averaging volume V which is centered at r. Note, that upon volume integration of the total density (14) the total dislocation line length contained within the volume is obtained. Analogously, integration of the signed curvature density (17) yields the number of closed loops as multiple of 2π. The vector of signed GND densities contains densities of positive and negative edge and screw dislocations as later used in the CDD evolution equation (23). We remark that the (scalar) geometrically necessary dislocation density ρ G is the norm of , i.e. ρ G ≡ = , and that the average line direction of the geometrically necessary dislocations is given by the unit vector l = / .

Numerical discretization of the averaging scheme
Numerically, the fields (14)- (17) can be discretized by replacing the spatial points r with a finite set of points r i which have the coordinates (x i , y i , z i ). On a regular grid each point r i is then the center of a sub-volume V i = x y z which defines the averaging sub-domain Ω i : and determines the resolution of the continuum representation. As a consequence of this coarse graining, it is well possible that inside a sub-domain Ω i line segments "cross each other" without intersections: although these lines are located on different physical glide planes they can not be resolved separately after coarse graining †. The discretized version of the scalar dislocation density (14) is: In (19) we compute the line length contribution of each dislocation segment by numerical quadrature: letting u k and w k be the abscissas and weights of the quadrature method of choice ‡ in the unit interval [0, 1], respectively, each quadrature point contributes a line length J(u k )w k to the domain Ω i which contains c(u k ). Conveniently, the quantity J(u k ) can be directly computed for each dislocation segment from (5) and (6). Analogously, † The coarse graining resolution is the reason why short-range interactions have to be handled differently as compared to DDD. This, however, will not be discussed here. ‡ In this work we have used a uniform distribution of the quadrature points i the interval [0, 1].
we define the discretized vector of signed GND densities, and the discretized signed curvature density as

Continuum Dislocation Dynamics
The Continuum Dislocation Dynamics theory provides conceptual steps to map a discrete dislocation system onto a set of continuous, density-like field variables together with the respective evolution equations governing their change in time. The original, higher-dimensional theory (hdCDD) was defined in a configuration space which added an additional dimension to the spatial domain [32,48]. hdCDD is very accurate and contains very detailed microstructure information [39,49] but suffers from requiring a large number of computational degrees of freedoms. As a computationally favorable albeit somewhat less accurate theory a simplified version (CDD) is used, which is based on spatial evolution equations for the total dislocation density ρ, a vector of geometrically necessary signed densities and the curvature density q [35,45,50]. We assume that the components of are oriented such that they represent the orientation of screw and edge dislocations: = [ s , e ]. The temporal evolution of these field quantities are -in local slip system coordinates -governed by the following set of equations: where n denotes the slip plane normal. The by 90 • rotated GND density vector is ⊥ = [ e , − s ], and we assume Q (1) = − ⊥ q/ρ (see [50] for a discussion of this assumption). The evolution equations are mathematically closed by an expression for the tensor A (2) which here is obtained from a "maximum entropy" approach [50]: where l ⊥ is the unit vector perpendicular to l . As shown by Monavari et al. [50] Φ can be approximated as Φ ≈ ( /ρ) 2 (1 + ( /ρ) 4 )/2 and is a non-linear interpolation between = 0 (isotropic dislocation arrangement) and = ρ (fully polarized dislocation arrangement). In section 6 we numerically investigate the validity of this approximation. The corresponding non-dimensionless scaling of the CDD equations (22)-(24) can be found in Appendix D of [51].
In general, the velocity v in the above equations needs to be specified in terms of a constitutive equation through which elastic interactions enter the model. This is a priori not part of CDD: e.g. long-range (or "internal") interaction stresses have to be obtained from the additional solution of a dislocation eigenstrain problem [42,52,53] while elastic short-range interactions have to be recovered in an alternative way based on CDD fields (see e.g. [54][55][56] for steps into this direction) due to the limited resolution of any continuous description. Within the present work we will not tackle this problem of "dynamic closure" and assume for benchmark purposes velocity fields given as stationary analytical functions that do not depend on the dislocation state §. In other words: we ignore all dislocation interactions and concentrate exclusively on the "kinematics", that is, how curved and connected lines move and evolve in space and time without interactions.

Regularization
Eqs. (19)-(21) can be used to extract coarse grained fields of geometrical data from a given configuration of discrete dislocations. However, the CDD evolution equations contain spatial derivatives which have to be approximated numerically and which thus require field data with a certain level of spatial smoothness. It is therefore helpful to replace the Dirac delta-function δ(r) in the convolution integral (7) with a regularization function G(r). To improve the numerical efficiency during computation of the convolution it is beneficial to make use of the sifting property of the Dirac function: we can always write the regularized function as G = G * δ, where the symbol * indicates convolution over three-dimensional space. Together with the fact that convolution and volume averaging commute, it follows that regularized fields can be obtained by convolution of (19)- (21) with G. In our numerical implementation, we chose a Gaussian standard distribution as the regularization function: where the standard deviation s characterizes the width of the dislocation density distribution. For numerical reasons G is approximated by a discrete Gauss function G d which is non-zero only at a finite number of discrete points r i and is normalized such that x y z i G d (i) = x,y G(r(x , y , z ))dx dy dz within a small numerical tolerance. We then compute the convolution of the coarse grained fields (eqns. (19)- (21)) with the discrete Gauss function G d . We note, that in (26) the only free parameter is the standard deviation s which subsequently will be chosen as the average dislocation spacings (the mean separation between a dislocation and its nearest neighbors) inside the averaging volume (see also the discussion in [53]). Note that choosing a much larger value for s would destroy details of the dislocation structure to a large extend while if s were very small this would result in a quasi-discrete rather than continuous density structure. This is further discussed and demonstrated in the appendix.

Numerical experiments
One of the key advantages of the CDD theory is the property of kinematic consistency, i.e. the ability to represent flowing dislocations through their characteristic geometrical properties: a dislocation is a linear, curved line where each of its points moves perpendicular to the respective local line tangents. In the following benchmark systems we will directly compare the evolution of discrete and continuous dislocation microstructure in analytically given velocity fields. We start with 3 situations of increasing complexity with regards to boundary conditions: dislocation loop distributions in a rectangular domain (i) with open boundaries allowing out-flux of dislocations, (ii) with impenetrable surfaces and (iii) with a circular obstacle (i.e. an internal boundary). All 3 systems will start from the same initial dislocation microstructure which makes it easy to compare the effects caused by different external and internal boundaries. The 4th benchmark system is a periodic system with a statistically homogeneous random distribution of dislocations and investigates how good fluctuation can be represented in a situation without the strong polarizing effect of boundaries.

System and initial values
The following examples study a system as sketched in Fig. 1 where a 3-dimensional pillar with a square base is discretized into thin sub-volumes along the z axis. We investigate one of such sub-volumes with size L x = L y = 2000b and height z = 75b, where b = 0.256 nm is the norm of the Burgers vector. Initial values for DDD simulations for Systems 1-3 are obtained by randomly distributing N d = 50 circular, discrete dislocation loops with radius R = 75b inside the volume V (note, that we only chose circular loops for ease of implementation; loops could be arbitrarily shaped -as long as they are closed). Each loop is positioned on a different glide plane and thus loops do not intersect by construction; for the first three systems their centers are distributed on a quadratic domain 0.5L x × 0.5L y such that all loops are completely contained inside the computational domain. The average density in this volume is ρ 0 = 4.8 × 10 15 m −2 . z is the coarse graining height such that the system can be numerically represented as a two-dimensional object -a thin lamella in the x − y plane (Fig. 1A). Continuous field data is obtained by 'smearing-out' the coarse grained fields through the regularization outlined in section 4.2. Initial values for System 4 are obtained by distributing the loop centers on the whole quadratic domain L x × L y while simultaneously taking care of the periodicity. For determining the standard deviation for the Gaussian distribution we estimate the mean dislocation spacing from the density ρ 0 , which results in a mean For our examples this estimate of the standard deviation works well for all investigated systems and at all time steps (further details on the choise of s can be found in the appendix). Only when the mean dislocation spacing changes drastically -in section 5.4 due to dislocation pile ups at grain boundaries -we adjust this value locally during the post-processing step for obtaining reference DDD field data to the actual value of the dislocation spacing. An example for DDD and resulting CDD initial structure that is used for all numerical examples is shown in Fig. 2: the total dislocation density has a maximum in the center region where most of the loops overlap. At the same time nearly equal numbers of loop segments of all possible orientations are present in the center region which makes the dislocations 'statistically stored'. Hence, the GND density | | tends to zero. The curvature density q has the same shape as the total density because initially all loops have the same radius and thus q = ρ/R. Although for each of the first 3 studies identical initial values are used, we note that the chosen CDD fields are not special or more suitable than any other set of statistically equivalent initial values.

Numerical methods
For the numerical solution of the continuum evolution equations a Galerkin finite element scheme together with an implicit time integration scheme is used. To accurately represent derivatives of up to second order, Lagrangian shape functions with cubic polynomials were chosen. The finite element mesh consists of 30x30 elements and thus has an element size which is in the same range as z. When high density gradients at boundaries occur (System 2 and 3) we additionally refine elements in those regions. To increase numerical stability in particular in regions with high fluxes and where the density approaches zero a small viscous term was added to each evolution equation. E.g. in (22) we replace the divergence term ∇ · (v ) with ∇ · (v + ∇ ), where is chosen sufficiently small such that the physical behavior of the equations is not affected while numerical oscillations are suppressed. The FEM simulations presented below all take less than two minutes computational time on a off-the-shelve computer with a 2.6 GHz dual core processor.
Computation of the evolution of plastic strain can be done for CDD based on the Orowan equation (∂γ t = ρbv) by time integration alongside with the CDD evolution equations (22)- (24). To obtain the initial values for γ from DDD one can compute γ from the GND vector and the relation = − 1 b ∇γ: geometrical construction of the swept area for each unit of GND density and subsequent superposition of these areasfollowed by division by z -gives the plastic strain. Since the plastic strain is a derived quantity it suffices to show that main CDD field quantities match those from DDD and we do not explicitly show γ.
The DDD model used the representation by chordal splines (α = 1); for further numerical details on the used Parametric Dislocation Dynamics method please refer to [14,46].

System 1: distribution of dislocation loops, open boundaries
For this system a constant velocity of v 0 = 0.01 nm/µs is prescribed, and time is expressed as multiple of the duration that density needs to completely traverse the domain, L x /v 0 . Boundaries are 'open', i.e. dislocations may leave the volume as if the domain were just a sub-domain of a much larger domain. This case of dislocation loops which have the same curvature and expand with constant velocity is a case for which the CDD theory is supposed to give mathematically exact results, and the system therefore can be used as a reference system for testing the quality of the numerical scheme and if the chosen standard deviation s is appropriate. Initial values are shown in Fig. 2. Two different simulations are run: discrete microstructure is evolved with the DDD model and -separately -the continuous microstructure is evolved with the CDD model. To compare the two simulations at time t the discrete microstructure was converted into CDD field quantities using the same strategy as for obtaining initial values from discrete data (in the sketch Fig. 1F which is then compared with E). Fig. 3 shows a snapshot in time, where in the contour plots the dashed lines indicate CDD values and the full lines are the contours of the converted DDD values. We observe that as dislocation loops expand they flow away from the center region where the total density ρ tends to zero (compare ρ in Fig. 3). The difference in the centers of the loops becomes more insignificant as their radii grow and the resulting density distribution takes a loop-like shape. In each point dislocations now have all approximately the same orientation, they become GNDs which shows in ρ ≈ | | as can be observed in Fig. 3. The comparison of DDD and CDD fields shows excellent agreement. Only minor features of the DDD fields are slightly smoothed out which is an artifact of the numerical solution. The volume integrals of ρ and q/2π give the total line length and the total number of dislocation loops and are shown in Fig. 4. For reference we also show the line length for an infinite size system (approximated by a significantly larger domain in x and y direction) for which the total line length can be obtained analytically as L(t) = N d 2π(R + vt). Initially, both systems exhibit a linear line length increase (loops expand with constant velocity). As soon as density leaves the volume through the boundaries the line length for the finite size system decreases; the line length in the infinite system, however, continues to increase linearly. CDD data matches the DDD data perfectly and also coincides for the large system with the analytical solution. The right plot of Fig. 4 shows the total number of dislocation loops. This value stays constant for the infinite system, which is correct since there are neither sources nor annihilation. For the finite system the surfaces act as sinks and the number of loops is reduced. Note that this interplay between density and curvature density will also become relevant when sources and annihilation are to be modeled (which we do not consider within the present work).

System 2: distribution of circular loops in a finite domain
Now the outflow of dislocations will be constrained by imposing impenetrable boundary conditions, i.e. the dislocation flux ρv is enforced to be zero at the boundaries. Numerically, this is done by prescribing a velocity field which decays smoothly to zero directly at the boundary (in physical terms this boundary layer could be imagined as the result of FIB damage in the surface-near regions). To mathematically model the boundary layer we use the following sigmoidal function, where the parameter β controls the shape of the function and is chosen such that the boundary layer is clearly visible and allows to study the dislocation flux in detail (β = 20). ξ is the normalized distance from the boundary, e.g. with ξ = 0 . . . 1. The velocity field can now be composed from the four contributions It is depicted in Fig. 5 and will be used for CDD as well as for DDD simulations. Both simulations are conducted with the initial dislocation structure that was used as well for System 1 (Fig. 2). Again, the DDD simulation starts from the discrete microstructure, the CDD simulation from the continuous microstructure as obtained by D2C conversion. Both simulations evolve independently. Only in a post-processing step we again convert the resulting DDD microstructure for reference purposes into continuous field quantities as outlined before. A snapshot in time of the respective  freely until parts of them reach the boundary layer within which the velocity decays to zero. The non-zero velocity gradient rotates line segments such that they eventually are parallel to the boundary, i.e. their line tangents are aligned with the contours of the velocity field. This causes the initial SSD density to increasingly become geometrically necessary when dislocations pile up against the boundaries (compare ρ and | | in Fig. 6). Additionally, dislocations need to bend with a high curvature radius near the corners of the domain. This shows in the high value of the curvature density in these regions. For a more detailed view Fig. 7 shows horizontal and diagonal line plots through the domain for two different time steps. We observe again an excellent agreement between the DDD and CDD field values for all times, and only small deviations can be seen close to the boundaries and in particular for the curvature density.  Fig. 6. The plotted data is along a horizontal line at y = L y /2 (blue) and a diagonal line at y = x (red), where the lower left corner of the domain is at x = y = 0 (cf. Fig. 1).
The total line length inside the system serves as a good plausibility test: initially we again would expect free loop expansion until finally all loops are aligned with the boundary and are (nearly) rectangular in shape. Fig. 8 shows the total line length vs time, and it is obvious that at all times DDD and CDD agree very well: both approach the horizontal tangent when all lines are aligned with the boundaries. The fact that DDD values are slightly below the horizontal tangent is correct because loops are not exactly rectangular but still have small rounded corners. For the DDD curve we directly used the line length from the discrete data while for the CDD data we first had to integrate ρ. The small difference in the final value of the line length arises from discretization errors due to very steep gradients at the boundaries. In this respect, more realistic systems that also consider elastic dislocation interactions will numerically behave even better because pile ups always have a finite size.

System 3: circular dislocation loops in a domain with an impenetrable inclusion
The third system studies the idealization of a plastically deforming matrix into which an elastic, circular inclusion is embedded. Dislocations can move freely in the matrix while the inclusion acts as an impenetrable obstacle for dislocation motion. Annihilation of lines is not considered, hence, the number of dislocation loops stays constant; for real systems this is a strong restriction, but again, the philosophy is to decouple interactions and reactions from the purely kinematic effects. The computational model is -as in System 1 -a quadratic domain (L x = L y ) with open boundaries and contains the inclusion in the center. The inclusion itself is represented by an internal boundary which we model as a circular sub-domain (radius R = L x /10) with zero velocity in the center region and a smooth transition from v = 0 to a constant velocity v 0 away from the inclusion (cf. Fig. 9(a)). As a result of the velocity field dislocations are first slowed down on their way towards the inclusion and eventually freeze in their motion so that they cannot enter the inclusion. Both simulation methods use the same initial values as before. This does place some dislocations inside the inclusion which is unphysical. It turned out, however, that the influence on the resulting microstructure is only small. Therefore, we decided to consistently use the same initial values as for all other systems as well. Evolving the DDD and CDD system gives typical dislocation patterns as shown in Fig. 10, additional plots including error plots can be found in the appendix.
An auspicious feature is the formation of pile-ups of bent lines around the inclusion: when expanding dislocation loops touch the inclusion they get a 'dent' (see the colored line in Fig. 10(a)) and adjust to the shape of the inclusion. At the same time, line segments further away start 'spiraling' around the inclusion and deposit an increasing number of segments around the inclusion (also see the figure in the appendix and the supplementary movie). This behavior has a number of consequences and can be observed in the CDD representation as well: • Circular pile-ups of dislocation density can be seen in the upper middle and right plots of Fig. 10 : ρ has a high value in the vicinity of the inclusion because density from piled-up dislocations together with pre-existing line segments in the inside of the precipitate add up. In conjunction with the GND density | | we infer that dislocations in the inside are mainly statistically stored while pile-ups around the inclusion are -literally -geometrically necessary.
• Circular GNDs around the inclusion form an inverted loop, cf. Fig. 10 bottom right and middle plot. Taking a look at the sign of e.g. the edge GNDs e around the inclusion it is found that their orientation was rotated by 180 • (also compare the initial values in Fig. 2). For instance the upper half is now negative and therefore would move downwards but is hindered by the inclusion.
• Bent dislocations around the inclusion have negative curvature -the above statement about an inverted loop already suggests that the curvature of dislocations around the inclusion should be negative which in fact can be observed in the two plots on the right of Fig. 11 where the average curvature k = q/ρ was computed and plotted along two lines passing through the inclusion. The inclusion has a radius of L x /10 which is equivalent to a curvature of k inc = 0.020 nm −1 . Due to the smooth interface the velocity is zero only at about L x /16 which is equivalent to a curvature of k min = 0.031 nm −1 . In the bottom right plot of Fig. 11 we find the negative curvature peak values at |k CDD | ≈ 0.040 nm −1 in good agreement with both the geometry value k min and the DDD curvature value of |k DDD | ≈ 0.028 nm −1 .
Deviations between DDD and CDD are larger than for the previous two benchmark models. They are particularly large in those regions which exhibit large gradients in the velocity and in the CDD field values, both of which cause increasing discretization errors. What is the average response of this system? Comparing the temporal evolution of the DDD line length in Fig. 10 to System 1 shows that the inclusion creates additional density which is clearly visible from approximately t = 0.15L x /v 0 on and is caused by lines bending around the inclusion. CDD is lagging behind and is initially even slightly below the line length of System 1 (due to the lower average velocity as compared to System 1). After t ≥ 0.25L x /v 0 the CDD model is then also showing an additional increase in line length.
The difference between CDD and DDD is the result of larger errors close to the  Fig. 10. The plotted data is along a horizontal line at y = L y /2 (blue) and a diagonal line at y = x (red), where the lower left corner of the domain is at x = y = 0 (cf. Fig. 1).
inclusion where CDD almost always underestimates the true density values. Increasing the resolution of the spatial discretization of the finite element scheme does not yield any appreciable alleviation. The reason is the very complex deformation state of dislocations that brings the CDD theory to its limits because near the inclusion the line curvature can be different for different line orientations in the same averaging volume -a detail that cannot be represented with this simplified variant of CDD (this CDD theory only can represent one average curvature value for each point/averaging volume). Although a number of details of this system have been predicted properly -e.g. the line curvature around the inclusion which is important when it comes to hardening effects in terms of line tension -this particular CDD formulation only roughly predicts the line length production close to the inclusion. More elaborate evolution equations which contain more information about the orientation distribution and curvature of dislocations are required as e.g. those derived in [45,50]. The D2C strategy for more detailed CDD formulations, however, will still remain the same.

System 4: random dislocation distribution in a periodic domain
The previous benchmark systems have in common that they start with random initial values which during time evolution become strongly polarized through the geometrical constraints imposed by boundary conditions (e.g. dislocations adjust their line orientation and curvature to the shape of external or internal boundaries and become geometrically necessary). System 4 will now investigate how good CDD performs in a bulk-like situation where the total density consists to a large extend of statistically stored dislocation density superimposed with (smaller) GND density fluctuations. Periodic boundary conditions are used to eliminate geometrical constraints. Initial values consist of a statistically homogeneous random distribution of 200 dislocation loops of the same radius, and their centers are distributed across the whole volume L x × L y × z. This results in the same average density ρ 0 as in the center region of Systems 1-3; all other parameters stay the same.  The velocity and loops' curvature radius are constant everywhere. Since this study analyzes the quality of the representation of flow and evolution of density fluctuations, it is useful to decompose the total density into fluctuating and volume-averaged contributions, ρ(x, y) = ρ fluct (x, y) + ρ . For plotting the field data we additionally normalize the fluctuations with the variance of the respective field, e.g. in Fig. 12 we plot for the total density (ρ − ρ )/var(ρ). Fig. 13 shows line plots of initial and evolved DDD and CDD data, additional error plots can be found in the appendix. Clearly, all details of the discrete microstructure are accurately reproduced. Fluctuations evolve in a non-trivial manner and we attribute negligible deviations to numerical errors. The average behavior matches exactly the DDD values. The overall very good agreement is not surprising after seeing that already System 1, which initially had a similar random structure in the center of the domain, gave very good agreement between DDD and CDD. It is, however, important to show that this good agreement is sustained during time evolution and that even the increase in the range of the fluctuations are properly represented. An application based on this type of system has been studied in the context of dislocation patterning [41] where in particular the fluctuations together with the correct line length increase and a Taylor-type flow stress were identified as key ingredients for pattern formation.

Determining the kinematic closure parameter Φ from DDD
The numerical studies in section 5 demonstrated that DDD simulations can serve as an elegant means for directly benchmarking continuum dislocation microstructure and its evolution. Thus, it seems natural to directly use DDD simulations to verify or even obtain expressions that were used to constitutively close the set of CDD evolution equations. Within the CDD equations occurred in particular one quantity which had to be assumed: the interpolation function Φ, which was taken in the 'maximum entropy approximation' approach as Φ = (| |/ρ) 2 (1 + (| |/ρ) 4 )/2. Φ was in an earlier approach assumed to interpolate linearly between SSD and GND density. We will now analyze DDD simulations in order to decide which assumption is under realistic simulation conditions the most suitable one. Recall the definition of the evolution equation of the curvature density and A (2) : In A (2) the approximation Φ occurs in a vector product together with a spatial velocity gradient. Thus, ∂ t q CD is sensitive w.r.t. variations in Φ only in regions where the velocity field is non-constant, as e.g. in the boundary layer close to the impenetrable boundary of System 2 from section 5.4. Continuous CDD field variables can be obtained from a DDD simulation as outlined before, and the time derivative ∂ t q can be numerically approximated by an explicit Euler scheme in time, At the same time also the evolution equation for q (eqn. (29)) can be evaluated based on field data obtained from DDD (Q (1) and ∇v are both known, and l, l ⊥ are functions of and ρ) . The unknown value of Φ can now be obtained for each point (where ∇v = 0) from solution of the inverse problem : Find Φ(x, y) -for all (x, y) with ∇v(x, y) = 0 -such that Note that no CDD simulation needs to be done, only DDD data is used. For each point (x, y) we record the optimum value of Φ together with the local total density ρ(x, y) and GND density (x, y). Plotting these data gives the distribution shown in Fig. 14. There, the point of time was chosen such that already a significant amount of density reached the boundary layer but did not yet reach a stationary state. For an earlier point of time dislocations would not yet have reached the region where ∇v = 0 and no data could be obtained. At a later point of time most dislocations form pileups of geometrically necessary configurations near the boundaries which results in a strong accumulation of data points at the top right corner of the diagram. It is obvious that the optimization strategy yielded points with a relatively large scatter. However, it can be concluded that most data points are located below the straight dashed line which is the linear interpolation between GNDs and SSDs. The polynomial maximum entropy approximation fits much better. Binning the data e.g. in 20 /ρ intervals gives averages that for /ρ = 0.8..1.0 coincide nicely with the polynomial approximation.
We note, that the space dependent tensors Q (1) and A (2) also can be extracted from DDD data. E.g. the '11' component of A (2) can be obtained by integrating the discrete density for each separate line segment against the cos 2 ϕ c of the discrete line orientation ϕ c , followed by regularization based on the Gaussian convolution.
Nonetheless, the data is not sufficient for any statistical analysis in all other regions and could not be much improved through additional simulations. These results might indicate that the approximation of ∂ t q is not sufficiently sensitive w.r.t. to the only parameter /ρ. Currently, higher order closure assumptions are under investigations. Their results might help to understand the difficulties in identifying a unique functional Φ for closure.

Conclusion
We developed a systematic method for averaging geometrical properties of discrete dislocation lines which allows for a direct comparison of discrete and continuum descriptions of evolving dislocation microstructures. Particular care was taken to formulate the numerical approximation such that discrete and continuous formulations are at all times consistent with each other (e.g. both the DDD and the CDD descriptions were based on the same spline representation). Using DDD data as reference the comparison of simulation results for 4 different systems revealed a surprising accuracy of the CDD theory and its numerical implementation: CDD is able to evolve very complex dislocation microstructure and simultaneously respect various physical boundary conditions in almost exactly the same level of detail as DDD simulations do. Even in situations where the simplifying assumptions of the used CDD formulation were clearly violated (System 3) we found that results were still reasonably. Within this work we only considered the kinematic closure and assumed a (fixed and analytical) velocity function. An important step for future work which still requires a large amount of fundamental work, however, is the "dynamical closure", i.e. the question what the average velocity of interacting continuum distributions of dislocations under stress is.
Considering that in 3D DDD simulations the computational cost scales approximately with N 2 (N being the number of nodes/segments) and considering that for a continuum theory the computational cost does not increase with the number of dislocations at all (not even when dislocation interactions are considered) we are very optimistic that a continuum description of dislocation plasticity might very soon complement DDD simulations in particular when it comes to high densities in large volumes and/or high accumulated plastic strains. Furthermore, the D2C conversion might even be directly beneficial for analyzing DDD data: the continuum representation is well suitable for ensemble averaging and could be a novel way of directly comparing discrete microstructures.

Appendix. On the resolution of continuum fields
In section 4 continuous fields were obtained that had to be artificially smoothed out for numerical reasons. This was done by convolving the coarse grained fields (eqns. (19)-(21)) with a discrete Gauss function G d . Choosing the standard deviation s in (26) is a choice which has to be made for every system and e.g. with regards to numerical aspects as well as based on the desired degree of detail which the density fields should be able to represent (e.g. in terms of density fluctuations). In the following, we summarize some key aspects that one might want to consider for deciding on a suitable value for s.
• The obvious criterion for a lower bound for s is related to the resolution of the numerical scheme: choosing s much smaller than roughly the size of a finite element (which additionally depends on the used shape functions) results in discretization errors, i.e. the fluctuations of the density field cannot be numerically represented properly anymore. In principle the discretization error can be measured. For practical purposes one could e.g. compare the numerical representation of a 'continuum loop' with the analytical solution. Note that choosing a value of the standard deviation in the range of a Burgers vector results in a computationally very expansive continuum model of nearly discrete objects (Fig. A1 left).
• A numerically stable solution of the transport equations of density fluxes additionally requires a certain degree of smoothness in order to avoid strong, undesired numerical oscillations. Studying the aspect of numerical stability is nontrivial (see e.g. [57] for a stability analysis for a continuum dislocation dynamics model based on the Kröner-Nye tensor) and was not attempted for the present CDD equations so far. We approach this problem in a pragmatic way and choose numerical parameters as well as the lower bound for the standard deviation s such that a single continuum dislocation loop can be expanded accurately within a small error tolerance during the simulated time [48].
• Deciding on an upper bound for the standard deviation s is problem specific: CDD works properly for all values of s exceeding the above introduced lower bounds. However, the larger the chosen s the more details of the heterogeneous microstructure will be destroyed. In particular, choosing s in the range of magnitude of the system dimension causes all gradients of the density fields to vanish. As a consequence, the partial differential equations for transport of densities deteriorate to ordinary differential equations, which are no longer able to describe fluxes on a scale below the system size (Fig. A1 right).
In this work, the emphasis is on heterogeneous dislocation microstructure. Hence, we chose, guided by the mean dislocation spacing of s ≈ L x /35, a value of s which is small enough that fluctuations are not getting smeared out, but which at the same time is large enough so that discrete dislocations can no longer be differentiated. Figure A1: Initial density ρ for different values of the standard deviation s obtained from the DDD data of System 4. The left plot is a nearly discrete microstructure, while in the right plot an almost homogeneous structure without any gradients resulted. CDD can in principle evolve all initial values accurately, the only limitation is our used numerical scheme, which does not work for values of below s ≈ L x /200.

Appendix. Additional time evolution and error plots for System 2, 3 and 4
In Fig. B2 -B4 the time evolution of discrete and continuous microstructure together with the relative error |ρ DD − ρ CD |/ρ DD for System 2 (grain with impenetrable boundaries), System 3 ('precipitate') and System 4 (periodic BCs) are shown. In order to reduce the unreasonable diverging behavior of the relative error in regions where ρ DD → 0 we chose a cut off value of 0.5 (i.e. less than half the size of the smallest interval of the color scale of the density contour plot for System 2 and 3) below which the error is not computed (shown as white regions). In the density plots the solid lines are -as before -the contours of the converted DDD data while the dashed lines show the CDD data. Fig. B2 has in most regions and for all time steps relative errors of less than ≈ 5%. Only at later time steps when dislocations pile up against the boundaries and density exhibits steep gradients deviations become larger in some places. The CDD theory does represent the evolution of this system properly. System 3 in Fig. B3 shows clearly that the geometrical constraint of the inclusion creates a complex situation which only partially can be represented by CDD and causes large errors. Despite the fact that a number of features are correctly represented (in particular the curvature close to the inclusion, see main text), a more precise evolution of this system would require a more refined CDD theory. System 4 in Fig. B4 demonstrates that CDD consistently and accurately predicts the increase of average density as well as the transport of density fluctuations -both of which are automatically contained within the set of evolution equations. The relative error is in most places even smaller than 4%. Figure B2: System 2 -DDD data, density and relative error at time steps t = 0, 0.2 and 0.4L x /v 0 (from top to bottom). The white regions in the error plots indicate that the relative error was not computed due to (nearly) zero density. Figure B3: System 3 -DDD data, density and relative error at time steps t = 0, 0.135 and 0.27L x /v 0 (from top to bottom). The white regions in the error plots indicate that the relative error was not computed due to (nearly) zero density.