Consistent determination of geometrically necessary dislocation density from simulations and experiments

The use of Nye's dislocation tensor for calculating the density of geometrically necessary dislocations (GND) is widely adopted in the study of plastically deformed materials. The curl operation involved in finding the Nye tensor, while conceptually straightforward has been marred with inconsistencies and several different definitions are in use. For the three most common definitions, we show that their consistent application leads to the same result. To eliminate frequently encountered confusion, a summary of expressions for Nye's tensor in terms of elastic and plastic deformation gradient, and for both small and large deformations, is presented. A further question when estimating GND density concerns the optimization technique used to solve the under-determined set of equations linking Nye's tensor and GND density. A systematic comparison of the densities obtained by two widely used techniques, L1 and L2 minimisation, shows that both methods yield remarkably similar total GND densities. Thus the mathematically simpler, L2, may be preferred over L1 except when information about the distribution of densities on specific slip systems is required. To illustrate this, we compare experimentally measured lattice distortions beneath nano-indents in pure tungsten, probed using 3D-resolved synchrotron X-ray micro-diffraction, with those predicted by 3D strain-gradient crystal plasticity finite element calculations. The results are in good agreement and show that the volumetric component of the elastic strain field has a surprisingly small effect on the determined Nye tensor. This is important for experimental techniques, such as micro-beam Laue measurements and HR-EBSD, where only the deviatoric strain component is measured.


Introduction
Before the development of crystal plasticity finite element formulations, phenomenological continuum models were used to describe plastic deformation in materials. Without delving into the underlying microstructural processes, these approaches, calibrated by experiments, could capture plasticity at the macroscopic scale (Dunne and Petrinic, 2004;Khan, Akhtar.S, 1995). Their great advantage is simplicity, however a lack of physical basis severely limits their predictive capabilities, especially for processes where microstructural heterogeneity is important.
Crystal plasticity finite element (CPFE) formulations address this issue by explicitly modelling plasticity in terms of crystallographic slip at the grain scale (Roters et al., 2010). Popularity of these formulations has increased dramatically as they directly account for complex interactions between individual grains of polycrystals and the resulting locally heterogeneous loading. Beginning from 1982, when it was first introduced by Peirce et al. (Peirce et al., 1982), the CPFE technique has developed to span a range of constitutive and numerical formulations, applicable to a large number of problems. For example CPFE has been used to simulate the development of microstructures and the consequent effect on the macroscopic material response (Aifantis, 1984), to simulate surface roughening in thin film mechanics problems (Raabe et al., 2003), grain-boundary and interface mechanics (Bate and Hutchinson, 2005;Meissonnier et al., 2001), strain-gradient effects (Dunne et al., 2012(Dunne et al., , 2007, polycrystalline morphology and texture, the necessary conditions (energy) for crack nucleation , geometrically necessary dislocations (GNDs) (Dahlberg et al., 2014), creep and high temperature deformation (Balasubramanian and Anand, 2002), texture formation (Asaro and Rice, 1977), deformation twinning (Kalidindi, 1998), multiphase mechanics (Vogler and Clayton, 2008) etcetera. Importantly CPFE methods can be easily adapted to different material systems, simply by modifying the crystallographic slip law. CPFE has been applied to a diverse range of not only metals (Balasubramanian and Anand, 2002;Dunne et al., 2012;Li et al., 2009;Vogler and Clayton, 2008), but also rocks (Behrmann, 1985). Furthermore CPFE has been used across a range of length scales, from single crystals (Wang et al., 2004;Weber et al., 2008) to polycrystals (Zhao et al., 2008) and multiscale applications combined with ab initio calculations (Raabe et al., 2007). CPFE has also been used with other modelling techniques such as continuum dislocation dynamics, and nonlinear thermoelasticity to simulate the response of materials under extreme dynamic loading (Luscher et al., 2016).
In the early 70s, empirical viscoplastic formulations were primarily used, based on the "flowpotential" approach proposed by Rice (Rice, 1971) for time-dependent plastic deformation. Subsequent works of Rice and Asaro (Asaro and Rice, 1977) and Peirce and Needleman (Peirce et al., 1982), focused on the analysis of non-uniform, localized deformation of ductile crystals, where crystal slip was simulated by a rate-independent, elastic-plastic relation, following the Schmid law.
Owing to computational restrictions, their simulations involved a simplified scenario of a single slip, or two symmetric slip systems.
With increasing computational power, a wide range of microstructure-based multi-scale plasticity models has emerged. These include various grain-and sub-grain scale problems, as well as complex 2D, 3D grain morphologies. The introduction of strain gradient terms in the constitutive law marked a major step forward, making it possible to capture experimentally observed size effects.
Numerous strain gradient plasticity (SGP) formulations have been proposed, for example by Fleck and Hutchinson (Fleck and Hutchinson, 1997;Fleck et al., 1994), Gao, Huang, Nix and Hutchinson (Gao et al., 1999), Arsenlis and Parks (Arsenlis and Parks, 1999), Cheong and Busso (Busso et al., 2000), Anand, 2005a, 2005b), Dunne et al. (Dunne et al., 2007) and Fleck and Willis (N A Fleck and Willis, 2009;N. A. Fleck and Willis, 2009). These approaches have enabled accurate simulations of inelastic, scale-dependent deformation phenomena such as the increasing hardness of metals and ceramics with decreasing indenter size in indentation simulations (Wang et al., 2004), or Hall-Petch grain size strengthening effect (Lim et al., 2014;Lyu et al., 2015;Raabe et al., 2003). SGP phenomenological formulations have also been extended to calculate the fraction of the rate of plastic work converted into heat, by taking into account a strain dependent factor to include the locked in strain energy around statistically stored dislocations (Lubarda, 2016).
In particular, SGP formulations helped to numerically simulate the length-scale effects and production of GNDs associated with non-uniform plastic deformation. Smaller characteristic lengthscales lead to steeper strain gradients and hence higher GND densities, causing a size effect as flow stress depends on dislocation density (Nye, 1953). Liu et al. (Liu et al., 2013) compared experimental and theoretical evaluations (using tension and torsion on polycrystalline copper wires) of three phenomenological theories of strain gradient plasticity, to show that the size effects seen in plastic flow is primarily due to the GND density generated as a result of plastic strain gradients. Using both mechanism-based and phenomenological SGP theories, Paneda et al. (Martínez-Pañeda and Niordson, 2016) showed localized strain hardening near crack tips, promoted by GNDs. Compelling as the results are, experimental techniques are required to confirm these numerical simulations. The critical thickness theory has recently been used to get a more reasonable estimate of the length-scale (in the µm range) from the underlying fundamental physical quantities to facilitate the use of the SGP theory in engineering applications such as finite element applications (Liu and Dunstan, 2017).
In strain gradient CPFE formulations, plastic deformation is accounted for by dislocation glide on active slip systems and the deformation gradient is linked to the lattice curvature and in turn to the additional GNDs, generated in the slip systems, to accommodate this lattice curvature. The lengthscale effect within the concept of GNDs is captured here by using Nye's formulation of the dislocation tensor (Nye, 1953). In a recent study, Lyu et al. (Lyu et al., 2015) modelled crystal plasticity using a continuum dislocation dynamics model (CDD) and used this together with a viscoplastic selfconsistent model to study the evolution of dislocation densities in multi-phase steels.
To calculate GND density, the closure failure of a suitable Burgers' circuit can be considered (Ashby, 1970;Nye, 1953). Using Stokes theorem, this can be recast as the computation of the "curl" of the deformation gradient to form the dislocation or Nye tensor (Nye, 1953). Whilst straightforward conceptually, this step has been marred with inconsistencies. A review of the literature reporting GND density calculations shows a wide range of different curl definitions being used, often with erroneous applications of ± signs, misplaced indices and missing transpose operations. These errors will lead to unphysical results. For example a missing transpose operation effectively corresponds to swapping of Burgers' vector and line direction, resulting in incorrect dislocation densities. If a minus sign is erroneously placed then left handed screw densities become right handed and edge dislocation densities have their extra half plane on the opposite side of the slip plane. The first key goal of this paper is to compare the three most commonly used curl definitions and to establish the correct expressions to be used.
A further question in the computation of GND density concerns the optimization technique used to solve the under-determined set of equations linking the curl of the deformation gradient and GND density. Two optimization techniques, L 1 and L 2 , are commonly employed. Each yields a different solution and very few studies (Wallis et al., 2017;Wilkinson and Randman, 2010) have investigated the differences in the results they produce. Here we carry out a systematic comparison of the GND densities predicted by both methods to determine their applicability in different scenarios.
For the validation of strain gradient CPFE models, a direct comparison to experiments performed at the same length-scale is essential. Here we present experimental measurements and strain gradient CPFE calculations of the lattice distortions beneath a spherical nano-indent in a tungsten single crystal. Experimental techniques such as electron back-scattered diffraction (EBSD), high-resolution EBSD (HR-EBSD), high-resolution digital image correlation (HR-DIC) are commonly used to measure lattice distortions in two-dimensions, for example studies by (Kysar et al., 2010), (Dahlberg et al., 2014), (Kartal et al., 2015), (Zhang et al., 2016), (Ruggles et al., 2016), (Guan et al., 2017) etcetera. Barabash et al. (Barabash et al., 2009) showed how GNDs and the effective strain gradient change the white beam Laue patterns of shocked materials. With the aim of capturing the GND formation in plastic deformation, we use the synchrotron X-ray micro-Laue diffraction technique to non-destructively probe the full lattice rotation and residual elastic strain field with 3D spatial resolution, and without altering the residual stress state. Using strain gradient CPFE calculations we carry out a detailed 3D simulation of the same experiment. The experimentally measured and predicted lattice rotations, strains and GND densities are compared in detail. This gives rise to several interesting questions, for example concerning the effect of the volumetric elastic strain on GND density calculations, since Laue diffraction only measures the deviatoric lattice strain tensor (Chung and Ice, 1999). This question is examined by modelling the experimentally recorded data using the strain gradient CPFE calculations.
We begin by describing the different definitions for calculating the curl of a second-order tensor. This is followed by a review of the theoretical framework of the computation of the dislocation tensor and GND density. Based on this the expressions for the dislocation tensor in terms of elastic or plastic deformation gradient, as well as lattice strains and rotations, are discussed. Next, a comparison of nano-indentation-induced lattice distortions measured by Laue diffraction and predicted by strain gradient CPFE simulations is presented. Finally, the effects of L 1 or L 2 optimisation techniques, and volumetric elastic strain on the computation of GND densities are explored.

Computing the Curl of a second-order tensor
As noted by Robert W. Soutas (Soutas-Little, 1999), there are several different definitions in use for computing the curl of a second order tensor. Here, three different approaches to the curl computation are discussed. Importantly we show that, if used consistently, they all lead to the same end result.
Let P and V be general second-order tensors. The km component of the pre-curl of is denoted as (∇ × ) , while the post-curl is ( × ∇) (Soutas-Little, 1999). In component form these may be stated as where, ∈ is the permutation tensor. An alternative third definition (referred to as "curl3" throughout this text) is commonly used in computational crystal plasticity studies, e.g. by Arsenlis and Parks (Arsenlis and Parks, 1999) and Cermelli and Gurtin (Cermelli and Gurtin, 2001): The derivation of each of these three curl formulae, in tensor notation, is shown in Appendix A.
The post-curl definition is the negative transpose of curl3. This can be shown as for any second-order tensor , the curl computation using each of the three discussed conventions, may be equated as

Calculation of dislocation tensor
A geometrical link between the lattice curvature and the distribution of GNDs is given by the dislocation tensor (also known as the Nye tensor), , proposed by Nye (Nye, 1953). Nye initially formulated using only the lattice rotation gradients, assuming that no long-range elastic strain fields are present. Kroner and Ashby further developed this approach, by adding elastic strain gradients to the formulation of the dislocation tensor (Arsenlis and Parks, 1999). Below a summary of this theory is provided.
The analysis below, for small deformations, closely follows the derivation by Fleck and Hutchinson (Fleck et al., 1994) in their study establishing the concept of strain gradient plasticity. Figure 1b shows a representative crystal lattice within an imaginary solid, with a chosen Cartesian reference frame as depicted. We assume that plastic flow occurs by dislocation motion and that the lattice is stretched and rotated during elastic deformation. Let us consider the relative displacement, , of two material points which are separated by as shown in Figure 1a. This relative displacement can be split into three components (Fleck et al., 1994) as follows = where, is the relative displacement caused by slip and described by the slip tensor . is caused by lattice rotation and is due to elastic strain, . For a specific slip systems, λ, defined by slip direction and slip plane normal , and crystallographic slip, , the slip tensor, , is given by the contribution from each of the active slip systems. Thus, Following Nye's reasoning (Nye, 1953) the closure failure of a Burgers' circuit, c, on surface S with plane normal N (see Figure E.1 (b)) can be used to link the crystallographic slip to the resultant Burgers' vector, < > : Using Stokes' theorem this can be rewritten as: where is the dislocation tensor defined by Nye. Thus, the Nye tensor corresponds to the curl of the slip tensor (Fleck et al., 1994): Nye (Nye, 1953) So, Nye's dislocation tensor may be written as where λ is a general slip system.
Here, it is important to note that the total displacement, , along any closed circuit must be zero.
Thus, the closure failure brought about by crystallographic slip, < > = ∮ , has to be balanced by an equal and opposite displacement incompatibility, i.e. ∮ ( ). This relation between the closure failure due to slip (i.e. plastic displacement) and that due to elastic displacement can be described by the concept of the deformation gradient as outline below: Figure 1a shows an imaginary material in the undeformed configuration with a line vector dX. After deformation, this line vector is transformed to dx. From here on we distinguish between undeformed and deformed coordinate systems. Variables in upper case correspond to undeformed coordinates, while lower case refers to the deformed coordinates.
The deformation gradient, F, can be defined as a second-order tensor that maps the undeformed state to the deformed state of a sample. This can be written as where, u is the displacement, is the identity matrix and the displacement gradient. contains information about whether the deformation involves a rigid body rotation or a shape change or both (Dunne and Petrinic, 2004). When the deformation includes both elastic and plastic contributions, the deformation gradient can be split into elastic ( ) and plastic ( ) parts, where each of the gradients may contain both stretch and rigid body rotation.
In Eq. (14) , is defined as a curl computation of the slip tensor. The plastic deformation gradient captures the deformation by crystallographic slip, which is the same as the deformation captured by the slip tensor (Eq. (8)). In fact, = . Hence Nye's dislocation tensor may be written in terms of as = ∈ , .
Kroner (Kroner, 1955) and Bilby (Lazar and Pellegrini, 2016) expressed the dislocation density tensor as the negative of the expression adopted by Nye (Eq. 14). Consequently, Arsenlis and Parks (Arsenlis and Parks, 1999) and Cermelli and Gurtin (Cermelli and Gurtin, 2001) used the curl3 convention to find the dislocation density tensor from the plastic deformation gradient. Rewriting (Eq. (14)) in terms of the curl3 formula gives From here onwards, unless otherwise specified, the curl operation signifies performing the curl computation using the curl3 convention. Owing to the contribution made by several researchers to the concept of the dislocation density tensor, from here on we adopt the notation of to represent it.
Eq. (18) can also be arrived at considering the following approach: To separate the elastic ( ) and the plastic ( ) deformations, an intermediate imaginary configuration dp can be introduced ( Figure   1b), where the sample has undergone purely plastic deformation. The transformation of dX to dp is captured by the plastic deformation gradient ( ) dp can then be mapped to the vector dx by the elastic deformation gradient = and equating Eq. (19) and (20) Rewriting Eq. (21) the multiplicative decomposition of the deformation gradient is obtained as proposed by Lee et al. (Lee, 1969) Substituting and into Eq. (17) and introducing elastic and plastic parts of the displacement gradient, and respectively, gives For small deformations, the term is negligible and thus the displacement gradient for small strains may be written as From the kinematics of deformation, closure failure of a region can be defined as the change in length of a path on the surface due to generation of dislocations in the volume. Acharya and Bassani (Acharya and Bassani, 2000) defined closure failure with respect to deformed configuration x, where where <b> is the net Burgers vector of the dislocation lines passing through closed loop c. This is analogous to Eq. (12) in the deformed configuration. Using Stokes theorem, the integration around c may be replaced by integration over any surface patch s, bounded by c and with plane normal n, so where, curl −1 ≅ −curl for small elastic strains because as curl of the identity matrix is zero. Re-writing in terms of the displacement gradient, The transpose in Eq. (28) is introduced when applying Stokes' theorem to higher order tensors, as proved by Cermelli and Gurtin (details in Appendix C) (Cermelli and Gurtin, 2001). The closure failure, represented in terms of the non-vanishing cumulative Burgers' vector of all dislocations, can also be written in terms of the undeformed configuration X. Computation of curl in the undeformed configuration, dX, will be referred to as "CURL" from here on.
Here <b> and <B> refers to the resultant Burgers' vector in the deformed and undeformed coordinate frame respectively. Rewriting Eq. (29) in terms of the plastic displacement gradient, Thus, in summary, the closure failure can be represented in terms of either the elastic or plastic displacement gradient as where for small deformation we do not need to distinguish between the initial and deformed system.
The elastic strain, , and lattice rotation, , are related to the displacement tensor by Thus, Eq. (32) may further be re-written as for the case of small deformation.
Given the deformation gradients, or , or displacement gradients, or , can be computed using any of the three curl definitions discussed above (Eq. (1), (2) & (3)). Following Eq. (6), a summary of expressions for in terms of the three curl definitions is provided in Table 1.
In terms of: In terms of: Table 1 -Summary of computation of dislocation tensor using the three common different definitions of curl in terms of and .
For small strains the dislocation tensor may be explicitly written in terms of the plastic ( ) or elastic ( ) displacement gradients: The displacement gradient can also be written in terms of the lattice rotation and lattice strain as This leads to the following form for the dislocation tensor: where, considering the asymmetric nature of , the components are 0. In general, the lattice rotation gradients are substantially greater than the elastic strain gradients and make a more significant contribution to the dislocation tensor.

Calculation of GND density
Eq. (16) where, is the line direction. Since generally j>9 there is no unique solution for . Instead, knowing and A, optimization methods may be used to obtain . The mathematically simplest is the L 2 optimization scheme (Arsenlis and Parks, 1999), which minimizes the sum of squares of dislocation densities i.e ∑ 2 = . . Using the right pseudo inverse, the solution may be written as When using the L 2 optimization, it is essential to construct including all possible slip systems not just the active ones (i.e. the calculation is independent of the resolved shear stress where, ν is the Poisson's ratio. The "linprog" algorithm, implemented in Matlab (The Mathworks Inc.; www.mathworks.com) was used to perform the L 1 optimisation.
A key assumption is that dislocations are either pure edge or pure screw. Other methods for solving Eq. 40 involve minimising the total dislocation density (Demir et al., 2009;El-Dasher et al., 2003;Sun et al., 2000), or minimisation the equivalent line length (Wilkinson and Randman, 2010).
A thorough, in-depth comparison all these minimisation norms is beyond the scope of this paper.
Instead we focus on a comparison of the two most commonly used methods, L 1 and L 2 minimisations.
For both these optimisations, the solution obtained represents only one of the infinite number of solutions to Eq. (39). The total dislocation density is obtained by summing the magnitudes of the densities of all j dislocation types.

Material and Methods
To illustrate the concepts discussed above, we now consider the lattice distortions, dislocation tensor and GND density beneath a spherical nano-indent in a pure tungsten single crystal. A direct comparison is made between experimental measurements and numerical predictions from a straingradient CPFE model of the indentation process.
A [001]-oriented high purity tungsten single crystal (99.99 wt.%) was mechanically polished using diamond paste and colloidal silica to produce a near defect-free mirror finish. 500 nm deep indents were made using a MTS NanoXp indenter with a spherical, ~4.2 μm radius diamond tip. Synchrotron X-ray micro-beam Laue diffraction was used to probe the residual lattice distortions beneath a specific indent with sub-micron (~0.5 microns) 3D resolution. Briefly, micro-beam Laue diffraction measurements were carried out at beamline 34-ID-E, Advanced Photon Source, Argonne National Lab, USA. A polychromatic X-ray beam (7-30 keV) was focused by KB mirrors to a probe spot of ~500 nm full width at half maximum at the sample. The sample was placed at this probe spot in 45° reflection geometry and the orientation of the laboratory coordinates in relation to the initial  single crystal tungsten cube (20 ×20×20 µ 3 ) representing one quarter of the experimental setup with elastic properties as stated in Table 2. Crystal plasticity was implemented in the model using a UMAT subroutine (details of the UMAT and the crystallographic slip law used is provided in Appendix E) and assumptions of isotropic elasticity and small deformations were made in the numerical simulation.  (Ayres et al., 1975;Bolef and De Klerk, 1962;Featherston and Neighbours, 1963;Klein and Cardinale, 1993).
The boundary conditions imposed on the tungsten block included symmetric boundary conditions on the XZ and YZ surfaces near the indent, a traction free top surface, and fixed displacement and rotation boundary conditions on the remaining surfaces. The modelled spherical indenter (4.2 µm radius) was assumed to be a discrete rigid part, and contact between the tungsten block and the indenter was defined using the Abaqus node to surface contact algorithm. Consistent with the nanoindentation experiment, in the simulation, a displacement of 0.5 µm was applied to the indenter. A refined finite element mesh (applied edge bias 0.1 to 2 µm) with >15700 20-noded, reduced integration (8 integration point) 3D quadratic elements was used (C3D20R). The experimentally measured nano-indentation load-displacement data was used to refine the critically resolved shear stress parameter (CRSS), used in the slip law, to ensure accurate reproduction of the loaddisplacement curve. The effective modulus, Eeff, (Table E.2) was taken into consideration (Li et al.,1 With the assumption of an isotropic, linear elastic solid, the Young's modulus and Poisson's ratio are related to the elastic constant as follows: = 11 − 2 ( 12 2 11 + 12 ) and = 12 ( 11 12 ) ⁄ . gradient crystal plasticity was implemented with a user material subroutine (UMAT) that shares data between gauss points using a common block. The UMAT code was based on the original user element developed by Dunne et al. (Dunne et al., 2007). Further details of the constitutive law and model are provided in Appendix E.
The model was constructed with the initial crystallographic orientation of the sample. Both, experimental and the simulated results are presented in the same sample coordinate frame to enable a direct comparison (Figure 4 and Figure 6).

Residual elastic lattice strains and rotations
Lattice orientation of all sample points, captured by rotation matrix R, was measured experimentally by Laue diffraction and also predicted by the CPFE simulations. The average of the rotation matrix, R, of points located between 22-25 µm beneath the indent (approximate location of the red dot in However, in their study, a two-dimensional deformation state was purposely introduced to eliminate all out-of-plane deformation gradients (within experimental error), such that the resultant dislocation tensor had only two non-zero components. Using the advanced technique of Laue diffraction, we have been able to measure the out-of-plane components of the lattice distortions and therefore the elastic portion of the deformation gradient with sub-micron resolution.
As seen in Figure 4, CPFE predictions match well with the experimental results, for all three rotation components, except at the indent centre ( Figure 4, slice 2), where a rapid variation of lattice rotations is seen. A quantitative comparison between the CPFE and the experimental results are made in Figure   5 where line plots corresponding to the contour plots in Figure 4    Same axis applies to all graphs Laue which is related to the total strain tensor ( ) by To make a direct comparison with the experimental measurements, the dilatational and plastic strains were removed from the total strain predictions from CPFE. Figure 6 shows the direct components of predicted by CPFE and as well as the experimentally measured strains, plotted on the same YZ sections at different positions along the X-axis. Qualitatively there is quite good agreement, especially for the ɛzz out-of-plane strain component. A quantitative comparison between the results is made through comparing line plots (Figure 7), corresponding to the contour plots in Figure 6 (b) and (d), extracted at depth of 5 µm beneath the indent (white dotted lines in Figure 6 (a) and (c)). The agreement of the lattice strains is not as clear as for the lattice rotations, but similar features can be identified in the measured and predicted profiles. This is especially the case for slices 1 and 3, i.e.
slightly (5 µm away) away from the indent centre. Figure 8 shows all the components of the symmetric deviatoric strain tensor plotted on XZ and XY planes through the indent centre. Apart from the shear components, ɛxz and ɛyz, the strains predicted by CPFE and those measured experimentally agree quite well. In particular the symmetry of the deformation fields is reproduced. Experimental data for the ɛxz and ɛyz components is noisy as the experimental configuration is relatively insensitive to these strain components (Hofmann et al., 2013). Same axis applies to all graphs Laue The key limitation of HR-EBSD and HR-DIC is their lack of depth-resolved information, thereby allowing examination of the deformation field only at the sample surface. In contrast, micro-beam Laue measurements allows 3Dresolved strain measurement with very good sensitivity of ~10 −4 .
However, its spatial resolution (~0.5 to 1 µm in 3D) is lower than HR-EBSD (~0.05 µm). Our measurements show surprisingly good agreement between the measured lattice rotations and strains, and those predicted by CPFE (Figure 4 and Figure 6), inspiring some confidence in the use of this combination of techniques for analysing crystal scale deformation. α was then found by taking the curl of .  ].

Dislocation tensor computation
An important question concerns the error incurred by neglecting the volumetric component of the elastic lattice strain. Consider the difference, δ, between (−curl ( )) and (−curl ( )) , which both provide approximations to : The difference, δ, corresponds to the curl of the volumetric strain component ( ): Figure 10 shows plots of δ, computed from the CPFE simulations, on the same YZ, XZ and the XY planes through the indent centre as used in Figure 9. A comparison of Figure 9 (column 3: depicting the components of calculated using (−curl ( )) ) and Figure 10 shows that the magnitude of δ is substantially smaller than that of . This suggests that effect of (−curl ( )) on the calculated components of is small. The 3D depth-resolved measurements of deviatoric lattice strain and rotation, possible with microbeam Laue diffraction, allow determination of all nine components of the dislocation tensor. This is in contrast to surface techniques, such as HR-EBSD (Wilkinson and Randman, 2010), (Wallis et al., 2016), (Jiang et al., 2015), (Ruggles et al., 2016) or micro-Laue diffraction without depth resolution (Irastorza-Landa et al., 2017), where terms of the dislocation tensor Eq. (38) depending on 3 remain unknown. Hence, without depth-resolution, only three of the nine elements of can be explicitly determined. If the assumption is made that the effect of lattice strains is negligible, five components may be determined (Pantleon, 2008). This means that GND densities determined from 2D surface methods will always constitute a lower bound estimate.

GND density computation
GND density was computed using L 2 and L 1 optimisation techniques for both experimental measurements and CPFE simulations. For tungsten we assume that deformation is accommodated by dislocations with a/2<111> Burgers' vector slipping on {110} planes (Marichal et al., 2013;Srivastava et al., 2013) (list of Burgers' vectors and line directions in Appendix G). Furthermore, we assume dislocations to have either pure edge or pure screw character. This results in 16 distinct dislocation types; four screw types with <111> line directions and twelve edge type with <112> line directions.

L 1 vs L 2 Optimisation
The GND densities of all sixteen dislocation types, determined using the L 2 optimisation method (Eq. (39) and (41)), are shown in Figure 11 and Figure 13 for experiments and CPFE simulations respectively. In both figures dislocation densities are plotted on YZ, XZ and XY sections through the indent centre. GND densities, determined using the L 1 optimisation method (Eq. (39) and (42)), are shown in Figure 12 and Figure 14 for experiments and CPFE simulations respectively. For the L 2 optimisation the GND density is distributed (almost evenly) over all slip systems. The reason is that minimisation of ∑ 2 = . associates a larger penalty with slip systems that have high dislocation density. Thus a solution where dislocation density is distributed amongst slip systems is favourable.
In contrast, the L 1 scheme minimises the total energy (weighted line length) namely Here a much greater variation of dislocation density distribution between slip systems is observed. Figure 11 -Experimental dislocation densities obtained by L 2 optimization method plotted at the indent centre on the YZ, XZ and the XY plane. Colour scale shows log10(ρ) with ρ in 1/µm 2 . Scale bar = 5 µm.   Arsenlis and Parks (Arsenlis and Parks, 1999) compared L 1 and L 2 optimisation techniques for an fcc crystal and found that L 1 produced more accurate results for a dislocation structure consisting of two dislocation lines. In contrast, the L 2 method predicted complex dislocation structures with multiple dislocation lines. Randman et al. (Wilkinson and Randman, 2010) and Ruggles et al. (Ruggles et al., 2016) also observed that L 1 minimisation generates an uneven distribution of GND density over individual slip systems. However, in their study, no corresponding comparison was made to densities obtained using the L 2 method. Our direct comparison of L 1 and L 2 methods, for both experimental and CPFE datasets, is consistent with these observations. It highlights that L 2 optimisation leads to an unphysical spreading of dislocation density over many slip systems, making the use of the L 1 method for accurate estimation of GND densities on individual slip systems essential. This effect will be of particular importance for crystal plasticity simulations where a distinction between the cutting density and mobile density associated with particular slip system is made (Roters et al., 2010).
The total dislocation density (i.e. summed over all slip systems) computed through both methods is remarkably similar. This is shown in Figure 15 where the total GND densities (calculated using L 1 and L 2 optimisations) from experiments and CPFE are plotted on YZ, XZ and XY sections through the indent centre. Thus, if only the total dislocation density in the sample is required, either of the optimisation techniques may be used. In this case, the L 2 optimisation, which is far more straightforward to implement and is computationally cheaper, is preferable. Figure 15 -Total dislocation density plotted on YZ, XZ and XY sections through the indent centre. The GND densities predicted by L 2 and L 1 optimisation based on CPFE data are shown in (a) and (b) respectively. L 2 and L 1 results from the experimental measurements are plotted in (c) and (d) respectively. The colour scale shows log10(ρ) with ρ in lines/µm 2 .

Conclusions
The dislocation tensor captures the dislocation population required to accommodate inhomogeneous plastic deformation. It is linked to the density of geometrically necessary dislocations (GNDs), and can be equated to lattice deformation in the form of lattice rotations and lattice strains. This provides a very useful relationship between GND density and the lattice curvature it accommodates. Here we have provided a comprehensive review of the theoretical background of the computation of Nye's dislocation tensor and the underlying GND density in plastically deformed materials. Comparing CPFE simulations and X-ray diffraction measurements of lattice distortions associated with spherical nano-indents in tungsten, a number of important conclusions can be reached:  The relationship between different curl definitions, used to compute the dislocation tensor, has been explored. Table 1 provides a summary of the dislocation tensor in terms of both elastic and plastic deformation gradients, for cases of small and large deformations.
Importantly the different curl definitions in use, if applied consistently, all lead to the same result.
 Lattice rotations and lattice strains beneath a spherical nano-indent in a tungsten single crystal are considered. CPFE and synchrotron X-ray micro-beam Laue measurements show good qualitative agreement, particularly for the 3D distribution of lattice rotations.
 The contribution of (−curl ( )) , the curl of the volumetric part of the elastic strain tensor, to the dislocation tensor is small. This is an important result since many experimental strain measurement techniques, such as white-beam Laue diffraction and HR-EBSD, can only measure the deviatoric lattice strain tensor.
 L 1 optimisation recovers a more heterogeneous distribution of GND density over individual slip systems. In contrast L 2 optimisation distributes GND density almost uniformly over all slip systems. Thus, if GND density on specific slip systems is required, the physically-based L 1 minimisation should be used.
 The total GND density determined by either L 1 or L 2 minimisation is remarkably similar. Thus the computationally simpler L 2 optimisation may be used if only total GND density is required.

Appendix A
Derivation of the three common definitions of the curl of a second order tensor.

Micro-beam Laue diffraction
In micro-beam Laue diffraction, the recorded images correspond to the sum of the intensity scattered by the entire volume illuminated by the incident beam. Thus, the depth along the incident beam from which a specific diffraction signal originated is unknown. Hence, if several grains are illuminated simultaneously, or if there are large lattice distortions, the Laue spots become broadened and difficult to interpret. At the 34-ID-E instrument this limitation can be overcome by carrying out depth-resolved Laue measurements using the Differential Aperture X-ray Microscopy (DAXM) technique. Here, a ~50 µm diameter wire is scanned in small steps between the detector and the diffracting sample. The depth vs intensity profile for each pixel on the detector is calculated by subtracting the diffraction images from the consecutive wire position increments and triangulating using the wire edge and the line of the beam. A detailed description of the DAXM technique and the 34-ID-E instrument is provided elsewhere (Hofmann et al., 2013;Liu et al., 2010Liu et al., , 2004.
Measurements were done to a depth of 20 µm beneath the sample surface. Laue diffraction patterns contained 30+ peaks and were indexed and fitted using the LaueGo software package (J.Z. Tischler: tischler@anl.gov) to extract both lattice orientation and the deviatoric elastic strain tensor at each measured point in 3D space. The measured strain and rotation gradients were then used to calculate the dislocation tensor and GND density.

Appendix E 3D CPFE model
In the 3D CPFE model, the mechanical response of the tungsten BCC crystal under indentation was predicted using a constitutive law incorporating crystallographic slip. A brief description of the constitutive law, originally developed by Dunne et al. (Dunne et al., 2012(Dunne et al., , 2007, is provided here. Recalling Eq. (22), it is known that the deformation gradient F, splits multiplicatively into its elastic and plastic parts. Plastic deformation occurs on a slip system, λ, when the resolved shear stress is greater than the critically resolved shear stress (CRSS). can be defined in terms of the crystallographic slip (relative displacement of two slip planes separated by a unit distance), slip direction s and slip plane normal n. The crystallographic plane normals and the slip directions are updated as the crystal lattice undergoes deformation. For a finite number of slip systems, is given by the sum of the contribution of the slip systems to the resultant slip The rate of change of is thus where ̇ is the crystallographic slip rate on slip system λ. The velocity gradient L is given by The velocity gradient can be split into symmetric and anti-symmetric components to give the rate of deformation D and the continuum spin W respectively. The total rate of deformation can be written as a sum of the elastic and plastic rates of deformation as is computed using Hooke's law, while is approximated by the symmetric part of . can be written as where, the higher products are ignored for small deformations.
When the UMAT is called by Abaqus, the UMAT is provided with the deformation gradient ( Where ′ is the stiffness matrix rotated into the sample coordinate system. The UMAT used here solves these equations implicitly, i.e. all quantities are written at the end of the time increment and the stress is forced to converge back onto the yield surface within a tolerance of 10 −12 MPa using the plastic corrector term ′∆ . The reduction of this stress residual ( = +∆ − ′∆ ) is done using the Newton-Raphson iterative method. Further details can be found in (Dunne et al., 2007).
The physically based slip law used here determines the slip rate on the slip system by considering the thermally activated process of movement of dislocations, overcoming pinning obstacles. For a slip system with average dislocation glide speed 〈 〉 Burgers' vector magnitude , with q dislocations per unit area h-L as shown in Figure E.1(a), the crystallographic slip rate maybe written using the Orowan equation where, is the density of mobile dislocations. Figure E.1 -(a) Schematic diagram of a set of slip planes, viewed edge on, each comprising of a random distribution of dislocations on one slip system (b) One of the slip planes from (a), viewed in cross-section, with slip plane normal N and slip direction S. C represents the closed circuit path around the slip plane used to determine the Burgers' vector of the black cutting dislocations. An edge dislocation (shown in red), pinned by cutting dislocations from another slip system, where l is the distance between the pinning points and d the distance jumped by the dislocation on overcoming the pinning dislocations.
The thermal activation process influences the glide velocity by enabling the pinned dislocations to overcome the energy barriers produced by obstacles. Figure E.1 (b) shows an edge dislocation, pinned by cutting dislocations from another slip system, where l is the distance between the pinning points and d the distance jumped by the dislocation on overcoming the pinning dislocations. The average glide velocity is given by where, Γ is the rate of escape of dislocations given by where, ν is the frequency of dislocation jumps, G is the Gibbs free energy, k the Boltzmann constant and T the temperature in Kelvin. In terms of Helmholtz free energy (∆ ) and applied stress field τ, the Gibbs free energy is where, is the activation volume. Substituting this expression of ∆ into Eq. (E.8) and taking into account the forward and backward activation, and the critically resolved shear stress , the slip rate for a slip system λ, maybe written as, ̇ = ( ) 2 exp (− ∆ ) sinh ( sgn( )(| | − ) ) (E.11) depends on the spacing between the pinning dislocations . A lengthscale dependent mechanical response arises as the strain rate and subsequent plastic deformation reduces as the GND density increases. We assume for simplicity, where, is a coefficient representing the probability of pinning.
The values of the properties in the constitutive law have been acquired from literature and are listed in Table E.1.

Appendix F Expected lattice rotations
As the block deforms due to indentation, it is expected that lattice rotations about the X-axis, will be positive and negative (right-handed rotation) on either side of the indent centre (depicted by the dotted blue line). Similarly, it can be expected that θy will be negative in the region labelled before indent centre and positive after the indent centre. Rotations about the Z-axis are anticipated to be small.