A unification of finite deformation $J_2$ Von-Mises plasticity and quantitative dislocation mechanics

We present a framework which unifies classical phenomenological $J_2$ and crystal plasticity theories with quantitative dislocation mechanics. The theory allows the computation of stress fields of arbitrary dislocation distributions and, coupled with minimally modified classical ($J_2$ and crystal plasticity) models for the plastic strain rate of statistical dislocations, results in a versatile model of finite deformation mesoscale plasticity. We demonstrate some capabilities of the framework by solving two outstanding challenge problems in mesoscale plasticity: 1) recover the experimentally observed power-law scaling of stress-strain behavior in constrained simple shear of thin metallic films inferred from micropillar experiments which all strain gradient plasticity models overestimate and fail to predict; 2) predict the finite deformation stress and energy density fields of a sequence of dislocation distributions representing a progressively dense dislocation wall in a finite body, as might arise in the process of polygonization when viewed macroscopically, with one consequence being the demonstration of the inapplicability of current mathematical results based on $\mathrm{\Gamma}$-convergence for this physically relevant situation. Our calculations in this case expose a possible `phase transition' - like behavior for further theoretical study. We also provide a quantitative solution to the fundamental question of the volume change induced by dislocations in a finite deformation theory, as well as show the massive non-uniqueness in the solution for the (inverse) deformation map of a body inherent in a model of finite strain dislocation mechanics, when approached as a problem in classical finite elasticity.


Introduction
It is by now an accepted fact that the plastic deformation of metallic materials is primarily an outcome of the motion of dislocation line defects and that the evolving distribution of these defects, i.e., microstructure, plays a pivotal role in determining the strength and mechanical properties of such materials. In particular, there appears to be scientific consensus that the accumulation of Ashby's [Ash70] 'Geometrically Necessary Dislocations' (GNDs) leads to the phenomenon of size-effect [Ash70,FMAH94,SWBM93,MC95] in micron-scale bodies, in that a stronger overall stress-strain response is observed with a decrease in size of the specimen. Dislocations also form intricate microstructures under the action of their mutual interactions and applied loads such as dislocation cells [MW76,MAH79,MHS81,HH00] and labyrinths [JW84], often with dipolar dislocation walls, and mosaics [TCDH95]. The presence of such dislocation microstructures, in particular their cell size and orientation, has a strong influence on the macroscopic response of materials [HH12,Ree06]. While such dislocation mediated features of mechanical behavior from small to large length scales abound, commensurate theory and computational tools for understanding and predictive design are lacking, especially at finite deformation.
Several candidates for continuum scale gradient plasticity models at finite deformation [AB00b, APBB04, TCAS04, EBG04, KT08, MRR06, NR04, NT05, LNN19, NT19, MRR06, PHG19], including inertia [KN19a], are available with the goal of accurately modeling elastic plastic response of materials. These models can predict length-scale effects and some even account for some version of dislocation transport. However, to our knowledge, there is no continuum formulation that takes into account the stress field of signed dislocation density and its transport at finite deformation. Current versions of Dislocation Dynamics models [DNVdG03,IRD15] accounting for some features of finite deformation, reviewed in [AA19], are also not capable of computing the finite deformation stress and energy density fields of dislocation distributions. Appendix B presents a brief review of the vast literature on strain gradient plasticity (SGP) theories, as well as some models for computing static finite deformation fields of dislocations within a couple stress theory where dislocations are modeled as singular force distributions. This paper reports on the further assessment of a modeling approach grounded in the idea of continuously distributed dislocations to understand plasticity in solids at finite deformation, as it arises from dislocation motion, interaction, and nucleation within the material. Based on partial differential equations (pde), the microscopic model -Field Dislocation Mechanics (FDM) [Ach01, Ach04] -takes into account the stress field and energy distribution of signed dislocation density along with its spatio-temporal evolution at finite strains. The model is capable of dealing with several features of defects in the crystal lattice at the atomic scale and, through a commonly used filtering approach in multiphase flows [Bab97], provides a formal pathway for posing a space-time averaged pde model at the meso and macro scales [AR06, Ach07, Ach11] termed as Mesoscale Field Dislocation Mechanics (MFDM). Built on rigorous kinematical and thermodynamical ideas, MFDM allows for a study of finite deformation mesoscale plasticity rooted firmly in quantitatively identifiable links to the mechanics of dislocations.
A finite element based parallel computational tool based on the MFDM framework was verified and validated in [AZA19] and assessed in [AA19]. Here, we use this computational framework to study the various problems outlined in the abstract of this paper. In addition, we also demonstrate strong Bauschinger effects in our extension of J 2 plasticity theory as an outcome of inhomogeneity induced by boundary constraints to plastic flow, without the introduction of any ad-hoc model beyond an isotropic model of work hardening. This paper is organized as follows: after introducing notation and terminology in the remainder of this Introduction, Sec. 2 presents the governing equations of finite deformation MFDM. The staggered computational algorithm for the quasi-static framework is discussed in Sec. 3. Section 4 presents the results of the three illustrative problems mentioned in the abstract. We end with some concluding remarks in Sec. 5.
Significant portions of the description of the theory in Sec. 2 are common with [AZA19], a paper developed concurrently with this work. We include this content here for the sake of being self-contained, and because the theory being discussed is quite recent in the literature and not commonly known.

Notation and terminology
Vectors and tensors are represented by bold face lower and upper-case letters, respectively. The action of a second order tensor A on a vector b is denoted by Ab. The inner product of two vectors is denoted by a · b and the inner product of two second order tensors is denoted by A : B. A superposed dot denotes a material time derivative. A rectangular Cartesian coordinate system is invoked for ambient space and all (vector) tensor components are expressed with respect to the basis of this coordinate system. (·) ,i denotes the partial derivative of the quantity (·) w.r.t. the x i coordinate direction of this coordinate system. e i denotes the unit vector in the x i direction. Einstein's summation convention is always implied unless mentioned otherwise. All indices span the range 1-3 unless stated otherwise. The condition that any quantity (scalar, vector, or tensor) a is defined to be b is indicated by the statement a := b (or b =: a). tr(A) and det(A) denote the trace and the determinant of the second order tensor A, respectively. The symbol |(·)| represents the magnitude of the quantity (·). The symbol a en in figures denotes a × 10 n .
The current configuration and its external boundary is denoted by Ω and ∂Ω, respectively. n denotes the unit outward normal field on ∂Ω. The symbols grad, div, and curl denote the gradient, divergence, and curl on the current configuration. For a second order tensor A, vectors v, a, and c, and a spatially constant vector field b, the operations of div, curl, and cross product of a tensor (×) with a vector are defined as follows: In rectangular Cartesian coordinates, these are denoted by where ε ijk are the components of the third order alternating tensor X. The corresponding operations on the reference configuration are denoted by the symbols Grad, Div, and Curl.
I is the second order Identity tensor whose components w.r.t. any orthonormal basis are denoted by δ ij . The vector X(AB) is defined by [X(AB)] i = ε ijk A jr B rk . The following list describes some of the mathematical symbols used in this paper. C : Fourth order elasticity tensor, assumed to be positive definite on the space of second order symmetric tensors E : Young's modulus µ: Shear modulus ν : Poisson's ratio C e : Elastic Right Cauchy-Green deformation tensor I 1 (C e ) : First invariant of C e φ : Elastic energy density of the material ρ : Mass density of the current configuration ρ * : Mass density of the pure, unstretched lattice (·) sym : Symmetric part of (·) m : Material rate sensitivitŷ γ 0 : Reference strain ratê γ : Magnitude of slipping rate due to statistical dislocations for the J 2 plasticity model γ k : Magnitude of slipping rate due to statistical dislocations on the k th slip system for the crystal plasticity model n sl : Number of slip systems τ k : Resolved shear stress on k th slip system sgn(τ k ) : Sign of the scalar τ k m k , n k : Slip direction and the slip plane normal for the k th slip system in the current configuration m k 0 , n k 0 : Slip direction and the slip plane normal for the k th slip system in the pure, unstretched lattice g 0 : Initial strength (Initial yield stress in shear) g s : Saturation strength g : Material strength Θ 0 : Stage 2 hardening rate k 0 : Material constant incorporating the effect of GNDs on evolution of material strength η: Non-dimensional material constant in empirical Taylor relationship (∼ 1 3 ) for macroscopic strength vs. dislocation density : Material constant, with dimensions of stress × length 2 , used in the description of dislocation core energy density b : Burgers vector magnitude of a full dislocation in the crystalline material h : Length of the smallest edge of an element in the finite element mesh under consideration

Theory
This section summarizes the governing equations, initial and boundary conditions, and constitutive assumptions of finite deformation Mesoscale Field Dislocation Mechanics (details in Appendix A), the model that is computationally implemented, verified, and assessed in [Aro19, AZA19, AA19] and further evaluated in this paper. These governing equations are: In (1), W is the inverse-elastic distortion tensor, χ is the incompatible part of W , f is the plastic position vector, gradf represents the compatible part of W , α is the dislocation density tensor, v represents the material velocity field, L = gradv is the velocity gradient, T is the (symmetric) Cauchy stress tensor, and V is the dislocation velocity field. The upshot of the development in Appendix A is that if L p = 0 then the system (1)

Boundary conditions
The α evolution equation (1a), the incompatibility equation for χ (1b), the f evolution equation (1c), and the equilibrium equation (1d) require specification of boundary conditions at all times.
The α evolution equation (1a) admits a 'convective' boundary condition of the form (α × V + L p ) × n = Φ where Φ is a second order tensor valued function of time and position on the boundary characterizing the flux of dislocations at the surface satisfying the constraint Φn = 0. The boundary condition is specified in one of the following two ways: • Constrained : Φ(x, t) = 0, for x on the boundary and for all times t. This makes the body plastically constrained at the boundary point x which means that dislocations cannot exit the body through that point, while only being allowed to move parallel to the boundary at x.
• Unconstrained : A less restrictive boundary condition whereL p × n (2) is specified on the boundary, along with the specification of dislocation flux α(V · n) on the inflow part of the boundary (where V ·n < 0). In addition to this, for non-zero l, specification of l 2γ sd (curlα × n) on the boundary is also required, whereγ sd is defined in (11).
The incompatibility equation (1b) admits a boundary condition of the form χn = 0 on the external boundary ∂Ω of the domain. This condition along with the system (1b) ensures vanishing of χ when the dislocation density field vanishes on the domain, α ≡ 0. The f evolution equation (1c) requires a Neumann boundary condition of the form (gradḟ )n = (α × V + L p −χ − χL)n on the external boundary of the domain. The equilibrium equation (1d) requires specification of standard displacement/velocity and/or statically admissible tractions on complementary parts of the boundary of the domain.

Initial conditions
The evolution equations for α and f , (1a) and (1c), require specification of initial conditions on the domain.
For the α equation, an initial condition of the form α(x, t = 0) = α 0 (x) is required. To obtain the initial condition on f , the problem can be more generally posed as follows: determine the f and T fields on a given configuration with a known dislocation density α. This problem is solved by solving for χ from the incompatibility equation and then f from the equilibrium equation as described by the system where t denotes a statically admissible, prescribed traction field on the boundary. This determination of χ, f , and T for a given dislocation density α and t on any known configuration will be referred to as the ECDD solve on that configuration. Hence, we do the ECDD solve on the 'as-received' configuration, i.e. the current configuration at t = 0, to determine the initial value of f which determines the stress field at t = 0. For the dynamic case, an initial condition on material velocity field v(x, t = 0) is required.
The model admits an arbitrary specification ofḟ at a point to uniquely evolve f from Eq. (1c) in time and we prescribe it to beḟ = 0.
2.3 Constitutive equations for T , L p , and V MFDM requires constitutive statements for the stress T , the plastic distortion rate L p , and the dislocation velocity V . The details of the thermodynamically consistent constitutive formulations are presented in [AA19, Sec. 3.1]. This constitutive structure is summarized below.
Core energy density Υ (α) := 1 2ρ * α : α Table 1: Constitutive choices for elastic energy density, Cauchy stress, and core energy density. Table 1 presents the Cauchy stress expressions for the Saint-Venant-Kirchhoff material . It also presents the assumed constitutive form of the mesoscopic core energy density (per unit mass) for the material. Tables 2 and 3 present the constitutive assumptions for L p and V , respectively, for the crystal and J 2 plasticity models. Table 4 presents the governing equation for the evolution of material strength g for the two models. The use ofγ sd in (12) stems from the fact that isotropic (or Taylor) hardening is used for the evolution of strength on every slip system with equal initial values, i.e., where the function h is defined in (12) and (5) is a simple modification of standard latent hardening phenomenology assumed in classical crystal plasticity (see, e.g., [PAN83]). In this paper, we focus on constitutive assumptions using the J 2 plasticity model while the crystal plasticity counterparts are presented for completeness. Also, isotropic hardening is not a necessary condition for the formulation.
All material parameters, except k 0 and l, are part of the constitutive structure of wellaccepted models of classical plasticity theory. Our model requires these two extra material parameters beyond the requirements of classical theory. The length-scale l (with physical dimensions of length) sets the scale for the mesoscopic core energy to be effective, and k 0 (non-dimensional) characterizes the hardening rate due to GNDs.
We mention here that l, introduced in (7) or (9) as a dimensional consequence of including the core energy, is not responsible for producing enhanced size effects and microstructure in MFDM. Rather, the 'smaller is harder' size effect decreases with increasing magnitude of l since its presence reduces the magnitude of the α field and, consequently, reduces hardening (12).
Crystal plasticity

Numerical algorithm
The finite element implementation of the system of equations given in (1) has been discussed in [AZA19] where detailed numerical algorithms, verification, and validation exercises are provided. Here, we briefly describe the general flow of the algorithm for the sake of being self-contained.
Modeling material behavior through the use of MFDM requires the concurrent solution to the coupled nonlinear system of pdes (1). To efficiently solve the system (1) using a staggered scheme in each time increment as shown in [RA05,RA06] for the small deformation case, it is augmented with the rate form of the equilibrium equation to obtain the material velocity field in the domain. In the absence of body forces and inertia, the rate form of the equilibrium equation is given as The material velocity, v, obtained by solving (13), is used to discretely update the geometry at each time increment.
The augmented system ((1) and (13)) is solved numerically using the following discretization methods: the Galerkin FEM for the equilibrium equation (1d), its rate form (13), and the evolution equation (1c) for the compatible part of the inverse elastic distortion; the Least-Suares FEM for the incompatibility equation (1b); and the Galerkin-Least-Squares FEM for the dislocation evolution equation (1a). We now summarize the algorithm for the quasistatic case.

Algorithm for quasistatic loading case
The augmented system of equations ((1) and (13)) is discretely evolved in time. An efficient time-stepping criteria and an intricate 'cut-back' algorithm are implemented to ensure stable and accurate evolution of state variables. The algorithm used to solve the augmented MFDM system for the quasistatic case is as follows: • Given the material parameters and initial conditions on α and the prescribed traction t (most often vanishing), ECDD is used to obtain f , χ, and T on the configuration at t = 0.
• Given the geometry and state variables (f t , χ t , α t , V n , L p , and g n ) at time t n , ∆t n := t n+1 − t n is then calculated based on the time stepping criteria explained in [AZA19, Sec. 4]. The following steps are used to evolve the system in a time increment [t n , t n+1 ]: 1. The rate form of the equilibrium equation (13) is solved on the configuration Ω n to obtain the material velocity field v n . This velocity field v n is used to (discretely) evolve the geometry to obtain the configuration at time Ω n+1 .
4. f n+1 is determined as follows: (a) f is evolved from (1c) to obtain f n+1 on the configuration Ω n+1 . 5. Once the state at time t n+1 is accepted after checking through the cut-back criterion [AZA19, Sec. 4], the above algorithm is repeated to obtain the new state at t n+2 .
The fields obtained from Steps 3. and 4a. above suffice to define an approximation for the stress field at t n+1 , using the hyperelastic constitutive equation. However, this may not satisfy (discrete) balance of forces at each time step (because the current geometry obtained from the rate-form of equilibrium and the stress field under discussion need not necessarily be 'consistent' with each other in the sense of discretely satisfying force balance on Ω n+1 ); therefore we periodically use the discrete equilibrium equation to correct for force balance (Step 4b. above) .

Results and discussion
We use the parallel computational framework of MFDM developed in [AZA19] to solve three fundamental problems in finite deformation dislocation mechanics and small-scale plasticity.
Before proceeding to discuss results, we mention some details pertaining to our calculations. For all the results presented in this work, the input flux α(V · n) and curlα × n are assumed to be 0 on the boundary. Also,L p is directly evaluated at the boundary to calcu-lateL p × n. All fields are interpolated using bilinear/trilinear elements in 2-d/3-d, unless otherwise stated. The Burgers vector content of an area patch A with normal n is given by b A = A αn dA, where α denotes the dislocation density field in the domain. When the dislocation distribution α is localized in a cylinder around a space curve as its axis and this cylinder threads the area patch A, then we denote its Burgers vector by b. It must be noted that b is independent of any area patch A for which the intersection of the 'core' cylinder and the patch is entirely contained within the patch, this being a consequence of the fact that div α = 0. We refer to the magnitude of the Burgers vector, |b|, as the strength of the dislocation. We define a measure of magnitude of the GND density as [AA19] All algorithms in this paper have been verified to reproduce classical plasticity solutions for imposed homogeneous deformation histories by comparison with solutions obtained by integrating the evolution equation (14) for the elastic distortion tensor F e to determine the Cauchy stress response for an imposed spatially homogeneous velocity gradient history, L: where L p is defined from Eq. (9) or (7) with l = 0, andg is given by (12) with k 0 = 0.
A typical schematic of the basic geometry used in most problems (further details are mentioned as required) is shown in Fig. 1. The averaged T 12 component of the stress tensor on the top surface is denoted by τ . It is calculated by summing the tangential components of the nodal reaction force on the top surface and then dividing by the current area (line length) of the surface.Γ represents the applied strain rate. At any time t, Γ denotes the engineering shear strain and is calculated asΓ t.

Plastic flow in thin confined layers: comparison with experiment
Recent experiments by Mu et al. [MHM14,MCM14,MZHM16] report results of micropillars subjected to axial compression that contain a thin layer of ductile copper sandwiched within a stiff and brittle ceramic pillar. Confined simple shear loading conditions are generated when the thin Cu layer is inclined at 45 • to the pillar axis [MZHM16]. The interfaces between the Cu layer and the ceramic bulk may be assumed to be plastically constrained as the ceramic is brittle. The experimental results demonstrate a power-law relationship with negative exponent between the inferred applied shear stress on the Cu layer versus its thickness. This exponent was found to be ∼ −0.2 for the as-deposited material [MZHM16, As noted in [MZHM16, pp. 5-6] and by Hutchinson [Hut19], and posed as a fundamental challenge to all higher-order strain gradient plasticity theories which predict a size effect due to interface constraints, such models inevitably predict a scaling exponent ≤ −1, representing a much stiffer response than observed in experiment.
The situation above has prompted a further formulation of strain gradient plasticity [DO19] introducing a new fitting parameter and adapted to this single experiment, as well as reformulations [KN19b] of existing SGP models wherein ad-hoc relaxation of boundary constraints dictated by the basic theory is suggested to accommodate the observed behavior with the justification that "there is a maximum value of the magnitude of the plastic strain gradient at the layer boundaries that can be supported before plastic straining starts at the boundary." Moreover, it is understood from the kinematics of dislocation motion that plastic shearing parallel to a boundary is not constrained by constrained plastic flow boundary conditions [GN05] [AR06, Sec. 2.2.1].
In contrast, here (Sec. 4.1) we use the finite deformation J 2 MFDM computational framework to model the constrained simple shear of a thin polycrystalline metallic strip without any adjustment to the basic model to accommodate this specific experiment. We focus on the result for the as-deposited material and correspondingly use a reasonable value of the initial yield stress to reflect the presence of an initial statistical dislocation density in the material.
The scaling exponent, denoted by β, is evaluated by assuming a power-law relationship where c is a scalar which is constant for a given strain and boundary conditions, and H denotes the film thickness. Moreover, since MFDM accounts for the finite deformation stress fields for the GNDs and their spatio-temporal evolution coupled to the underlying kinematics, we also present the microstructure evolution during the deformation.
.3647 Table 6: Parameter values used to model the effect of thickness on the flow stress of thin metal films. Simulations are performed on rectangular domains of width 5 µm and varying thicknesses, H, of 0.5 µm, 0.65 µm, 0.8 µm, and 1 µm (the pillar diameter in the experiments was 5 µm). Problems are set up in a 2-d plane strain setting as follows: velocity boundary conditions corresponding to the simple shear deformation for plane strain condition are imposed. At any point x ≡ (x 1 , x 2 ) on the boundary denoted by P , v 1 =Γ y(x 2 ) and v 2 = 0 are imposed where y(x 2 ) is the height of the point P from the bottom surface as shown in Fig. 1. Table  6 presents the values of the material constants used in this section. The applied strain rateΓ = .001s −1 . The top and bottom surfaces are treated as (plastically) constrained whereas the other two lateral surfaces are treated as (plastically) unconstrained, as defined in Sec. 2.1. These boundary conditions mimic the metallic sample sandwiched between two non-deforming brittle substrates while undergoing shear deformation, and the possibility of free exit of dislocations from the sides of the strip sandwiched within the micropillar. Section 4.1.1 explores the effect when all the four surfaces are treated as plastically constrained. Figure 2 shows the convergence of the stress strain response for the domains with thickness 0.5µm and 1µm for two mesh sizes. We use the coarser mesh for our numerical computations, with details of the meshes for all the domain sizes given in Table 5. Figure 3 shows the scaling exponent β at different applied shear strains. The value of β predicted by MFDM is between 0 and −1 and very close to the value observed in experiments on the as-deposited samples in [MZHM16]. We predict a clear decreasing trend in the values of |β| with increasing strain. We note that the predicted values of β are a pure outcome of the theory unlike [DO19] wherein the scaling exponent is exactly equal to the fractional order of the discrete plastic-strain derivatives embedded in their theory, but undetermined by it. Figure 4 shows the ρ g distribution at different strains for the sample with H = 0.5 µm. In the figure, the x 2 axis is scaled by the thickness H of the sample. The constrained boundary condition on the top and bottom surfaces induce gradients in L p near these boundaries. The gradients, in the x 2 direction, of L p 21 -arising from the presence of the T 21 stress component -leads to the continuous generation of α 23 near the top and bottom boundaries. Figure 5 shows the ρ g distribution at Γ = 0.40 for the various sample sizes. The x 2 axis is scaled by the respective thickness of the sample. The figure shows that the non-dimensional dislocation layer width decreases with the increase in the height of the sample, an effect not predictable by classical theory on dimensional grounds [RA06].
In closing, we note that MFDM has not been shown to predict a size-effect at initial yield in small deformation analysis (at finite deformation, there is a variation of initial yield stress w.r.t the initial GND distribution field to be theoretically expected that needs to be assessed in future simulation work). Even though in the modeling of [DO19,KN19b] it is implicitly assumed that the entire size effect is related to phenomena at initial yield, it is not entirely clear to us from the experimental load-displacement data in axial compression in [MZHM16, Fig. 4b] that the plateaus there, representing axial displacements up to shearoff of the sample, do not include hardening in the stress-strain response. The results we have presented are directly related to the size-effect scaling in work hardening response. Our computational infrastructure is ideally suited for independent prediction of the finitedeformation load-displacement response of the experiment, which is future work in progress.

Effect of fully constrained boundaries and Bauschinger effect
Simulations are performed on square domains of sizes (1 µm) 2 , (2 µm) 2 , (5 µm) 2 , (10 µm) 2 , and (100 µm) 2 with a typical schematic of the geometry shown in Fig. 1. The boundary conditions are as for the simple shear case just discussed in this Sec. 4.1 except that the lateral boundaries are also considered to be plastically constrained. Table 6 presents the values of the material parameters used in this section. An applied strain rate ofΓ = .001s −1 is used. Figure 6a shows the convergence of stress-strain curves for the (1 µm) 2 sample size for two meshes with 70 × 70 and 140 × 140 elements. We use the coarser mesh of 70 × 70 elements for all the domain sizes. The unconstrained case represents a conservative simulation scenario with smaller gradients than the constrained case, and hence these mesh sizes suffice for the unconstrained case as well. Figure 6b shows the overall stress-strain response for all the domain sizes and both the boundary conditions (plastically constrained and unconstrained), demonstrating the 'smaller is harder' size effect.   Figure 7 presents the scaling exponent β at different strains. The values of β lie between 0 and −1 with magnitudes less than the case studied in Figure 3. The data over the whole size range does not fit a single power law expression; as shown, two power-laws appear to provide a reasonable fit. Figure 8a shows the stress strain plot for the (1 µm) 2 sample size with constrained boundaries up to 60% strain and Figure 8b shows the ρ g distribution at that strain. The loading direction is then reversed and the body starts unloading elastically. Fig. 8c shows the ρ g distribution in the domain at 59.17% strain when the averaged load on the top surface (τ ) is close to 0. As the reverse loading continues, the body starts deforming plastically again around 58.65% strain displaying a strong Bauschinger effect, presumably due to the internal stresses of the α distribution. Figure 8d show the ρ g distribution at 58.65% strain when the plastic deformation initiates again. Figure 9 shows the ρ g distribution in the same problem at three different strains of 40%, 45.46%, and 49.99% during the forward loading. There is considerable development of walllike microstructure in the interior of the domain between 40%-50% strain. This interior

Finite elastic fields in polygonization
We use finite deformation FDM to study the stress and energy density fields of a sequence of dislocation distributions whose limit is a through-dislocation wall, as observed in the physical process of polygonization [Gil55,Nab67]. After presenting our results, we make contact with available mathematical work [MSZ15b, Gin19b] on the limit energy functionals for nonlinear elastic deformations with dislocations, show that the current state-of-the-art of mathematical work in this direction is inadequate for the problem of polygonization, and discuss our perspective on the problem.
We compute the stress field of a special sequence of dislocations, wherein the dislocation cores stack up to form a dislocation wall in the limit. The result is a tilt grain boundary consistent with a piecewise uniform finite rotation field resulting in zero-stresses in the do- main. Before analyzing the stress fields for the sequence of dislocation distributions, we first study the limiting case, i.e., a polygonized domain with two dislocation walls as shown in Fig. 10. The height of the dislocation walls is chosen as 25k where k denotes the width of the walls.
x 1 In the limiting case, the orientation of the lattice on one side of the wall differs from the other side by a rotation angle denoted by θ 0 . To construct the dislocation density field describing the two walls in the domain for a given θ 0 , we first approximate this piecewise constant rotation angle field by a continuous field with the walls centered at x 1 = ±25k, and a is a dimensionless scalar chosen to be 0.238 which ensures that the widths of the dislocation walls are k. Figure 11 shows the variation of θ(x) along x 1 in the polygonized domain. The corresponding elastic distortion F e is then given by the rotation tensor field with the corresponding dislocation density field in the domain given by α = −curl W .
A non-uniform mesh, highly refined close to the dislocation walls, is used to discretize the polygonized domain comprising an isotropic Saint-Venant-Kirchhoff elastic material with E = 200GPa and ν = 0.3. The stress fields are calculated for θ 0 = 45 • in a 2-d plane strain setting. Since θ does not vary in the x 2 direction Thus, α 23 ≈ 0 is a reasonable approximation for low-angle grain boundaries. The solution for the nonlinear elastic stress field of the chosen F e field is of course trivial -it is 0 everywhere. However, the stress field of a progressively forming dislocation wall whose elastic distortion approaches a rotation tensor field everywhere is non-trivial, as we show in the following.
Since the latter is our main interest, and the knowledge of the trivial solution does not help in solving the latter problem, we devise a unified strategy to solve the entire class of problems. This requires first to solve the 'full-wall' problem with the trivial solution by a non-trivial method that then finds crucial use in solving the full range of 'sparse' to 'full' wall problems. In what follows, we first describe how the full-wall problem is dealt with.
With the dislocation density field and vanishing applied tractions specified, the ECDD system (2)-(3) is solved in the polygonized domain to determine the stress field. The numerical solution of the nonlinear ECDD problem requires a good initial guess. This guess is obtained by solving for f by using the Least-Squares finite element method with objective functional 1 2 Ω ||gradf + χ − W || 2 dV with W and χ specified. The specified data is generated from W , the inverse elastic distortion field corresponding to the prescribed F e field given by (17), and by solving for χ from the div-curl system (1b) with α = −curl W specified. The weak form (to solve for f ) is given by (18) f is fixed at one point in the domain to obtain a unique solution. Solving for f from the least squares method implies (W −χ) n = (gradf )n holds weakly on the external boundary. Figure 13 shows the distribution of the T 12 stress component in the domain, with the dislocation density fields shown in Figure 12. As can be seen, the shear stresses are negligible in the entire domain and further decrease upon refinement. To verify that the values of all the components of the stress tensor T are negligible, we define two non-dimensional measures of strain energy density in the domain as follows: In the above,ψ f d andψ sd denote the finite and small deformation non-dimensional strain energy densities, respectively. Figure 14 shows the distribution of finite deformation strain energy densityψ f d in the domain. The negligible magnitude of theψ f d distribution in   the body demonstrates that the body is stress-free. However, Figure 15 demonstrates that the small deformation theory predicts a non-zero strain energy density profile even when the values of the elastic distortion field is a rotation tensor everywhere. Moreover, the expression for the strain energy density for the small deformation (linear) theory is not invariant under superposed rigid body motions.
We now demonstrate the stress and energy density field paths induced by a sequence of dislocation distributions, going from one core to a full wall, in a rectangular domain of size 100k × 50k. The θ(x) corresponding to a dislocation wall centered at x 1 = 0 is given by where θ 0 is the rotation angle, equal to the difference in the orientations of the lattice across the wall. Equation (17) gives the elastic distortion tensor field in the domain. The values of a and θ 0 are chosen to be same as in the previous section, i.e. 0.238 and 45 • , respectively. As before, the negative curl of the inverse elastic distortion field gives the dislocation density comprising the full dislocation wall. The dislocation density field corresponding to a single core amounts to isolating the dislocation density distribution in a region of dimension k × k from the full wall, with value 0 everywhere else in the domain. The dislocation density corresponding to multiple cores in a wall configuration can then be prescribed by using the dislocation density field corresponding to the single core with a vertical shift and aligning the so-obtained core with the previous ones. For any given number of dislocation cores in the rectangular body, denoted by n c , the ECDD system is solved to obtain the stress field in the body. The chosen sequence of positions and number (of cores) for the α 13 component of the dislocation density is shown in Figure 16. α 23 is similarly placed and its values are accordingly assigned. The Burgers vector for a single dislocation core is computed to be b = .778k e 1 − .322k e 2 . A uniform mesh of (0.1k) 2 sized elements is used to discretize the rectangular domain for all cases except for the case of n c = 50. For n c = 50, the mesh is further refined near the dislocation wall.
For the simulations with n c ≤ 30 the initial guess to the Newton-Raphson method is obtained from the small deformation equilibrium problem solved on the current configuration as discussed in [ZAP18,Aro19]; in solving problems of finite deformation dislocation fields with sparsely distributed individual dislocations, this is found to be essential. However, for the simulations with n c ≥ 36, the initial guess is obtained by solving (18) (corresponding to a full dislocation wall), and this is also found to be essential to obtain a solution for the stress fields of these 'dense' distributions of individual dislocations. For n c = 35, the numerical solution did not converge with initial guess coming from either of the two approaches mentioned above. For n c = 32, a solution, using the initial guess from the small-deformation theory, can be obtained. The same procedure does not succeed for n c = 33, 34.    cores n c ≥ 36. The variation (for both T 12 and ψ f d ) which was spread out in the entire domain suddenly becomes localized near the wall and the boundary when n c ≥ 36. It is clear that when the number of (same sign) dislocations in the body increases from 0, the stress and energy have to increase. It is also clear that the full wall must represent a stress and energy-free configuration. Whether the transition between these behaviors has to be an abrupt 'phase transition' (w.r.t. the number of cores), as indicated by our calculations, is an interesting question for further study. When θ 0 is reduced, the magnitudes for the stress and energy fields become less pronounced with the qualitative conclusions remaining the same. Figure 19 shows the image of the constant e 1 vector field under mapping by the elastic distortion field on the current configuration for the full dislocation wall. We can see that the lattice on the left of the dislocation wall is rigidly rotated w.r.t. the lattice on the right by the prescribed misorientation angle θ 0 . We also calculate the change in volume (per unit length of the domain in the e 3 direction, i.e. the change in area) of the body due to the presence of the dislocation wall as described in Sec. 4.3. The % change in volume evaluates to 0 up to machine precision, a necessary condition for the elastic distortion field to be a rotation tensor field (which is spatially inhomogeneous in this instance).

Contact with the mathematical literature
In closing this section, we make contact with the mathematical results of [MSZ15b,Gin19b]. A summary of the main result of [MSZ15b] adequate for the present context is as follows: Given a frame-indifferent elastic energy density function W and ε a nondimensional measure of the magnitude of the Burgers vector of a dislocation proportional to a lattice constant of a crystal, 1 where F e ε is a sequence of elastic distortions consistent with a corresponding sequence of discrete dislocation distributions in the body given by µ ε through curl F e ε = µ ε , µ ε → µ, there exists a sequence of spatially constant rotation fields R e ε → R e such that R eT ε F e ε − I → β with |R eT ε F e ε − I| = O(ε| log ε|), and curl β = R eT µ, C is the standard linear elastic moduli (with minor (and major) symmetries), ϕ is a function that is deduced, and all Figure 20: Physical length scales for a single dislocation in a finite body arrows represent appropriate (scaled) convergences as ε → 0. Furthermore, an important hypothesis, for this discussion, on the 'well-separated'ness of admissible discrete dislocation distributions is that for any two dislocations comprising µ ε located at x and y, |x−y| > 2ρ ε , where "ρ ε ε, with ρ ε → 0 as ε → 0." That the sequence of rotation tensors R e ε with limit R e exists for the sequence F e ε (consistent with admissible µ ε ) is predicated on the assumption that the energy of the sequence as ε → 0 is bounded, i.e., Clearly, for the case of polygonization involving a single dislocation wall as we have considered, a representative sequence F e ε would have to converge to a rotation field that is not spatially uniform. Consequently, such a sequence would not be within the purview of the result of [MSZ15b], since the β for this sequence is neither skew-symmetric nor 0 so that the first term on the RHS of (22) does not vanish whereas the limit of the LHS of (22) does for this sequence. Thus, this failure has to be related to the assumptions behind the [MSZ15b] analysis. Since such hypotheses, and assumptions of the same ilk, are commonplace in the mathematical literature on dislocations [GLP10, Gin19a, Gin19b, LL16] and has begun to find prominence in the mechanics/engineering literature as well [RSC16,RDOC18], we discuss them here and provide our perspective on the matter.
First, we clarify an issue that is very rarely addressed in the literature (cf. [RSC16,RDOC18]), which, nevertheless, is crucial for understanding the physical setting of the mathematical results. This relates to the use of a fixed, bounded configuration containing dislocations whose strengths are assumed to scale with ε → 0. With reference to Fig. 20, we first note that the total energy of the cylinder of length L with a single screw dislocation is given by E = The energy content of the body can be interpreted in two different non-dimensional limits as ε → 0, corresponding to the expressions E µ2πLR 2 = ε 2 | log ε| and E µ2πLb 2 = | log ε|.
The first considers the body to be of fixed radius R with the interatomic spacing of the crystal b → 0 in which case the limiting total energy of a dislocation in the body vanishes.
With this physical interpretation of ε → 0, the ε 2 | log ε| 2 scaling employed by [MSZ15b] may be interpreted, but not necessarily, to correspond to a 'weakly unbounded' population of dislocations in the body growing as | log ε| as ε → 0 [GLP10, Gin19a, Gin19b] (we discuss a different physical interpretation of this 'fixed body scaling' in the next paragraph). We note, however, that for a physically realistic, non-degenerate scenario where a single dislocation in a body does not result in 0 internal stress and elastic strain energy regardless of its size, the interatomic spacing of a crystal is a finite, well-defined physical length for a specific material, as are the dimensions of a body, so choosing ε = b R → 0 with R fixed does not appear to be a viable proposition to us (when considering a fixed material).
On the other hand, keeping b fixed and sending R → ∞ seems eminently reasonable (the limit of a progressively large body of the same material), in which case the total energy of the body containing the single dislocation diverges as | log ε| → ∞ as ε → 0. We believe that it is this physically realistic scaling that should be employed in the analysis of dislocations and their distributions, with the minimum energy scale being sup ε→0 1 b 2 | log ε| Ωε W ( F e ε ) dV < ∞ (the domain depends on ε as it grows in size scaled by bε −1 and we refer to fields on the growing domains with an overhead ). By the nondimensionalization given by x = xε b with the corresponding domain of fixed size represented by Ω, it can be seen that It is in this sense that we interpret all 'fixed domain' dislocation-related asymptotic results with discrete dislocation strengths tending to 0, noting, in particular, the large difference between Ω W (F e ε ) dV and the physical energy content Ωε W ( F e ε ) d V as ε → 0. Consider now a 'full' dislocation wall described by a piecewise constant rotation field (as constructed in our computational example, but now with a sharp boundary) in a sequence of square domains of size H × H, with H steadily increasing, representing a progressively large domain containing a bicrystal. The magnitude of the Burgers vector of each individual dislocation in the wall is assumed to be b. Then the result of the case 1 ε = H b 1 may be approximated by the results of H b → ∞, b fixed. At least within the mathematical context being considered, the entire sequence of domains has 0 elastic energy (i.e. the part coming from W ), which is certainly O of any of the energy bounds ε 2 | log ε| or ε 2 | log ε| 2 as ε → 0. Moreover, it can be checked that, on the non-dimensional domain of fixed size, the strength of the Burgers vector of individual dislocations indeed scales like εb, if b is the magnitude of the Burgers vector of the dislocations on the physical sequence of domains. Why then does the elastic part of the limit energy in (22) fail to predict the energy content of such a wall, a commonly observed dislocation configuration, well worthy of prediction? The answer must lie in the fact that the 'full wall,' while being a very low-energy configuration, cannot be achieved as a limit configuration of the admissible dislocation sequences allowed by the [MSZ15b] analysis due to the 'well-separated'ness assumption -indeed, it is, in a sense, an 'opposite' limit that is considered in our computational example where the core dimensions 1 It is presumably based on this correspondence between the non-dimensionalized problem on the fixed domain and the growing domain problem that [MSZ15a,Gin20] state that "From the point of view of physics it is more natural to fix the lattice spacing and to consider domains 1 ε Ω of increasing size. Upon elasticity scaling both points of view are equivalent and fixing Ω rather than the lattice spacing is more convenient for the analysis. Thus ε really is a dimensionless parameter of the order of lattice spacing divided by the macroscopic dimension of the body. For brevity we will nonetheless often refer to ε as the lattice spacing." remains fixed and the inter-core distances reduce to 0. An interesting feature of this specific example is that the number of dislocations in a configuration (when the core radius tends to 0) does not correlate, even roughly, with the total energy content of the configuration, as might be expected based on energy scales related only to self-energy of single dislocations, an argument valid in relatively 'dilute' limits [GLP10,Gin19a,Gin19b]. Indeed, the defining energetic feature of the 'full wall' is that the entire energy of the individual dislocations in the wall gets screened by the '(nonquadratic) interaction energy' between them resulting in a very low-energy state.
Unlike Equation (2.4), Proposition 4.3, and Theorem 4.6 of [MSZ15b], the Generalized Rigidity Estimate (GRE) [MSZ15b, Theorem 3.3] does apply to the 'full wall' configuration under discussion. It is instructive to understand the breakdown of that result in validating the elastic part of the limit functional in (22) for the specific example of the full dislocation wall. As mentioned earlier in conjunction with (22), for admissible sequences (µ ε , F e ε ) satisfying the energy bound (23), |R eT ε F e ε −I| is small so that a quadratic approximation of W estimates it well and it is that quadratic approximation of W with argument β that appears in (22). For the 'wall sequence' (on the growing domains, and recalling that b H =: ε) with constant unit tangent to the wall in the the direction e 2 , | U e ε − I| is small (actually 0), where U e ε is the right stretch tensor of the polar decomposition of F e ε and the question becomes as to whether the spatially non-uniform rotation tensor field of F e ε can be approximated well by a constant rotation R e ε in the domain Ω ε so that U e ε can be well approximated by R e ε T F e ε . We are unable to deduce the necessary control on the point-wise values of R e ε − F e ε from the GRE, given here, for each fixed ε, by since, although the F e ε is a rotation tensor field for this specific sequence taking exactly two distinct values, say R e ε1 = R e ε2 , the various ingredients of the GRE in this specific example are given by and the constant C(Ω ε ) depends on the domain; for the corresponding problem on the fixeddomain, the point-wise squared magnitude of the difference R e ε − F e ε is therefore bounded at most by an O(1) quantity independent of ε (it is physically obvious that no theorem can prove that the rotation field of a high-angle, symmetric tilt boundary can be well-approximated by a constant rotation everywhere). In case such an approximation were to actually fail, then a 'quadratization' of W (about 0) is not a good approximation for it, and the limit functional in (22) cannot be the correct one for this wall sequence. Predictions of such a model would be similar in spirit to what is shown, e.g., in Fig. 15. Interestingly, however, the limit elastic functional in (22) is invariant under superposed rigid deformations.
The above arguments also seem to suggest that for small enough energy scales a better target for the elastic part of the limit functional is 1 2 Ω (U e − I) : C(U e − I) dV , where U e ε → U e as ε → 0, which is invariant under superposed rigid deformation as well as succeeds for the case of the dislocation wall. However, it should be noted that in the presence of large 'infinities' of dislocations of sufficient strength, control on the energy (and hence the elastic strain field) does not control the rotation field, this being expected since such control is a hallmark of compatibility of deformations, as in (non)linear elasticity theory without line defects.
Finally, in the context of energy scales of dislocation configurations for asymptotic analysis, the highest energy scales that have been considered, to our knowledge, are and sup ε→0 Ω in both cases, no limit energy functional is deduced as in the other references mentioned, and no (approximate) methods for computing stress fields of the dislocation distributions are devised.

Volume change due to dislocations
A fundamental question that was asked by Toupin and Rivlin [TR60] concerns the change in volume of a body when dislocations are introduced; linear elastic theory is not capable of capturing this volume change, due to the fact that the average value of each of the infinitesimal strain components of a self-equilibrated stress field in a body vanishes. However, it has been observed experimentally that the volume of the body changes upon the introduction of dislocations [Zen42] and the prediction by linear elastic theory (of no volume change) does not seem to be in agreement with experimental observations. Toupin and Rivlin [TR60] used a second-order approximation of nonlinear elasticity to give explicit expressions for the changes in average dimensions of elastic bodies resulting from the introduction of dislocations.
Here, we use finite deformation FDM to capture this volume change. The problem is set up in a 2-d plane strain setting as follows: Edge dislocations are assumed to be present in a rectangular body of dimensions [−10b, 10b] × [−10b, 10b]. An edge dislocation, with core centered at point p = (p 1 , p 2 ), is modeled by prescribing a dislocation density of the form The constant φ 0 is evaluated by making the Burgers vector of the dislocation core equal to be 1 , i.e. Ω α 13 da = b. The external boundaries are assumed to be traction-free. The body The volume change is calculated as follows: The ECDD system (2)-(3) is solved on the current configuration to obtain the inverse elastic distortion field for a given dislocation density in the domain. The % volume change is then calculated from (25) where V ref and V curr denote the volume of the reference and the current configurations, respectively: In the following, we present the volume change (per unit length along the e 3 axis) approximation for two cases corresponding to the total Burgers vector of all the dislocations being a) positive and b) zero. Table 7a shows the % change in the volume of the body consequent upon the introduction of multiple dislocations of the same sign, distributed in the body as shown in Figures 21a -21e. Table 7b shows the % change in the volume of the body due to the introduction of multiple pairs of dislocations, distributed in the body as shown in Figures  21f -21i, such that the total Burgers vector of all dislocations is 0. The configurations for zero resultant Burgers vector utilize positive and negative straight edge dislocation.
Tables 7a and 7b show that the change in volume upon introduction of dislocations as quantified by finite deformation FDM is very small for both the cases. Moreover, we see that this volume change is not linear w.r.t the number of cores i.e., the change in volume in the presence of 6 dislocations is not the same as 3 times the change in volume when 2 dislocations are present. The nonlinearity of the ECDD system and interaction among the dislocations cause the deviation from linear response to give smaller values.   This study suggests that for an isotropic Saint-Venant-Kirchhoff material, the volume change calculated with nonlinear theory is very small and therefore the linear theory prediction of no volume change is a very good approximation.

Non-uniqueness of inverse deformation in classical finite elasticity with dislocations
Consider the case when a dislocation with Burgers vector b is assumed to be present in the body. The solution to the ECDD system (2)-(3) on the current configuration Ω gives the inverse elastic distortion field W in the domain. It can be shown that if the dislocation core is removed from the body and a cut is made that runs from the boundary of the core to the external boundary of Ω to produce a simply connected body Ω s , then there exists a deformation y of this Ω s such that where W −1 =: F e and grad s denotes the gradient on Ω s . The reference configuration is then obtained by mapping this 'hollowed and cut' configuration Ω s by the field y. Moreover, the field y has two important properties: i) the jump in the value of y (denoted by y ) along the cut surface is equal to the Burgers vector b of the embedded dislocation in the original body Ω and, ii) y is independent of the cut surface chosen as well [ZA18].
Section 4.4.1 shows that the above topological properties are preserved in the framework of (computational) FDM. More interestingly, we are easily able to do similar calculations for a body containing multiple dislocations when the current configuration cannot be rendered simply-connected by a single cut, as shown in Sec. 4.4.2 -in this case we show that the jump in y is not constant along the boundary, corresponding to the cuts, of the simply-connected domain, in contrast to the single dislocation, single-cut case.

Inverse deformation for a single dislocation
We obtain non-unique reference configurations of a body Ω with a single dislocation with Burgers vector b = be 1 . The body Ω is made simply connected by removing the core and making a straight cut, at some angle from the e 1 axis, that runs from the boundary of the core to the external boundary. The field y on this simply connected configuration Ω s is calculated by using the Least-Squares FEM. Figure 22 shows the reference configurations for the cut surface chosen at varying angles w.r.t. the e 1 axis.
The reference configurations self-penetrate when the cut makes an angle between 0 • and 180 • with the e 1 axis. When the angle lies between 180 • and 360 • , the reference configurations show detachment. Also, the reference configurations shown in Figures 22b and 22c, and 22a and 22e are entirely different, even though the corresponding simply connected configurations for these cases are related by a rigid rotation of −/ + 90 • about the e 3 axis. This is a consequence of the constraint that the (vectorial) jump in y has to remain constant along the cut surface. Moreover, this jump is exactly equal to the Burgers vector be 1 of the dislocation in the original configuration Ω. The jump in the value of the field y is identical for all the cuts showing that the jump is independent of the cut surface, thus verifying the expected topological property.

Inverse deformations for multiple dislocations
This section explores the possibility of obtaining non-unique reference configurations of a body with multiple dislocations. Figure 23 shows such a scenario when 4 dislocations are assumed to be present in the body and reference configurations are obtained by solving for y using the Least-Squares FEM for different simply connected configurations. On the top of Figure 23, we show different simply connected configurations, each obtained by making multiple cuts (shown in red) in Ω. The corresponding reference configurations are shown on the bottom in Figure 23. As can be seen from these figures, the reference configuration is non-unique -for a given self-equilibrated W -as it depends on the cuts made in the body to make it simply connected. Moreover, even when the cut-configurations differ by rigid rotations (for example 23a and 23d, and 23b and 23c), the corresponding reference configurations (23e and 23h, and 23f and 23g,) are entirely different. Another important thing to note here is that unlike the case of the single dislocation, the jump in y is not constant along the cut in the presence of multiple dislocations, but is constant along the cut between any two cores, in accord with the analytical values they should take in these examples of configurations with multiple cuts.
The above examples highlight an important feature of the ECDD formulation -if the dislocation stress problem for specified dislocation density is formulated classically, i.e., involving only an inverse deformation field y on the 'deformed' configuration with the dislocation considered given, it is clear that there is massive non-uniqueness of solutions y for the same stress field corresponding to W 3 ; ECDD suffers from no such ambiguity and lends itself to robust computation.

Conclusion
A partial differential equation based model of finite deformation plasticity coupled to dislocation mechanics has been presented, with no restriction on material and geometric nonlinearities. The model fundamentally accounts for the stress field of arbitrary dislocation distributions and melds elastic dislocation theory at finite deformation with J 2 plasticity in a practical manner suitable for application. This paper along with [AA19, AZA19] form a trio that solves some key problems of current and classical significance in the fields of plasticity and dislocation mechanics, utilizing a finite element based computational framework developed in [AZA19]. As well, Lauteri and Luckhaus [LL16] state "It is worthwhile to compare our result with the differential geometric description of dislocation structures, introduced by Kondo, Kröner and Bilby at al. (see [Kon64,Krö59,BBS55]) and also with Γ-limit results in the context of linear elasticity where implicitly or explicitly a volume density of dislocations is assumed (see [GLP10]). It remains to be investigated if these models remain valid as averaged limits if on an intermediate scale there exists a Cosserat structure of micrograins." -up to the intended meaning of "averaged limits" above, we believe that we have answered this latter question in the affirmative, albeit with a far-reaching extension of the pioneering differential geometric works, going beyond simply kinematic considerations at finite deformations. We believe that our work provides a practical pathway for exploring many problems at the intersection of finite deformation plasticity and dislocation mechanics at realistic length and time scales. Further work will involve probing problems of shear localization in 2 and 3 dimensional bodies of rate-dependent and rate-independent materials, among others. Finally, improving the constitutive assumptions inherent in L p based on averaging of dislocation dynamics at the individual defect scale remains a fundamental pursuit. and since the theory being discussed is quite recent.
Field Dislocation Mechanics (FDM) was developed in [Ach01, Ach03, Ach04, Ach07] building on the pioneering works of Kröner [Krö81], Willis [Wil67], Mura [Mur63], and Fox [Fox66]. The theory utilizes a tensorial description of dislocation density [Nye53,BBS55], which is related to special gradients of the (inverse) elastic distortion field. The governing equations of FDM at finite deformation are presented below: Here, F e is the elastic distortion tensor, χ is the incompatible part of W , f is the plastic position vector [AR06], gradf represents the compatible part of W , α is the dislocation density tensor, v represents the material velocity field, L = gradv is the velocity gradient, and T is the (symmetric) Cauchy stress tensor. The dislocation velocity, V , at any point is the instantaneous velocity of the dislocation complex at that point relative to the material; at the microscopic scale, the dislocation complex at most points consists of single segment with well-defined line direction and Burgers vector. At the same scale, the mathematical model assigns a single velocity to a dislocation junction, allowing for a systematic definition of a thermodynamic driving force on a dislocation complex that consistently reduces to well-accepted notions when the complex is a single segment, and which does not preclude dissociation of a junction on evolution.
The statement of dislocation density evolution (27a) is derived from the fact that the rate of change of Burgers vector content of any arbitrary area patch has to be equal to the flux of dislocation lines into the area patch carrying with them their corresponding Burgers vectors. Equation (27b) is the fundamental statement of elastic incompatibility and relates the dislocation density field to the incompatible part of the inverse elastic distortion field W . It can be derived by considering the closure failure of the image of any closed loop in the current configuration on mapping by W . Equation (27c) gives the evolution equation for the compatible part of the inverse elastic distortion field. It can be shown to be related to the permanent deformation that arises due to dislocation motion [Ach04]. The field gradf can also be viewed as the gradient of the inverse deformation for purely elastic deformations. Equation (27d) is the balance of linear momentum (in the absence of body forces). Balance of mass is assumed to hold in standard form, and balance of angular momentum is satisfied by adopting a symmetric stress tensor.
The equations of FDM (27) can also be succinctly reformulated aṡ but since the system of Hamilton-Jacobi equations in (28) 1 is somewhat daunting, we work with (27) instead, using a Stokes-Helmholtz decomposition of the field W and the evolution equation for α in the form of a conservation law. Equation (27) is augmented with constitutive equations for the dislocation velocity V and the stress T in terms of W and α [Ach04, ZAWB15] to obtain a closed system. FDM is a model for the representation of dislocation mechanics at a scale where individual dislocations are resolved. In order to develop a model of plasticity that is applicable to mesoscopic scales, a space-time averaging technique utilized in the study of multiphase flows (see e.g. [Bab97]) is applied to microscopic FDM [AR06,Ach11]. For any microscopic field m given as a function of space and time, the weighted, space-time running average field m is given as where Ω is the body and Λ is a sufficiently large interval of time. B(x) is a bounded region within the body around the point x with linear dimension of the spatial resolution of the model to be developed, and I(t) is a bounded interval contained in Λ. The weighting function w is non-dimensional and assumed to be smooth in the variables x, x , t, t . For fixed x and t, w is only non-zero in B(x) × I(t) when viewed as a function of x and t .
Mesoscale Field Dislocation Mechanics (MFDM) is obtained by applying the above spacetime averaging filter to the FDM equations (27) with the assumption that all averages of products are equal to the product of averages except for α × V . The governing equations of MFDM [AR06, Ach11, AA19] at finite deformation (without body forces) are written as where L p is defined as The barred quantities in (29) are simply the weighted, space-time, running averages of their corresponding microscopic fields (defined in (27)). The field α is the Excess Dislocation Density (ED). The microscopic density of Statistical Dislocations (SD) at any point is defined as the difference between the microscopic dislocation density α and its averaged field α: which implies with b the magnitude of the Burgers vector of a dislocation in the material, ρ t the total dislocation density, ρ g the magnitude of ED (commonly referred to as the geometrically necessary dislocation density), and ρ s is, up to a scaling constant, the root-mean-squared SD.
We refer to ρ s as the scalar statistical dislocation density (ssd ). It is important to note that spatially unresolved dislocation loops below the scale of resolution of the averaged model do not contribute to the ED (α) on space time averaging of the microscopic dislocation density, due to sign cancellation. Thus, the magnitude of the ED is an inadequate approximation of the total dislocation density. Similarly, a consideration of 'symmetric' expansion of unresolved dislocation loops shows that the plastic strain rate produced by SD, L p (30), is not accounted for in α×V , and thus the latter is not a good approximation of the total averaged plastic strain rate α × V .
In MFDM, closure assumptions are made for the field L p and the evolution of ρ s , as is standard in most, if not all, averaged versions of nonlinear microscopic models, whether of real-space or kinetic theory type. As such, these closure assumptions can be improved as necessary (and increasingly larger systems of such a hierarchy of nonlinear pde can be formally written down for MFDM). In this paper, we adopt simple and familiar closure statements from (almost) classical crystal and J 2 plasticity theories and present the finite element formulation for the model. Following the works of Kocks, Mecking, and co-workers [MK81,EM84] we describe the evolution of ρ s through a statement, instead, of evolution of material strength g described by (12); L p is defined by (7) (or (9)) following standard assumptions of crystal/J 2 plasticity theory and thermodynamics. A significant part of the tensorial structure of (7) and (9) can be justified by elementary averaging considerations of dislocation motion on a family of slip planes under the action of their Peach-Köehler driving force [AC12].
Below, and in system (1) as well as the rest of the paper, we drop the overhead bars for convenience in referring to averaged quantities.
As shown in [AZ15], (1a) and (1b) implẏ up to the gradient of a vector field, which is re-written as where F e := W −1 . This can be interpreted as the decomposition of the velocity gradient into an elastic part, given byḞ e F e−1 , and a plastic part given by F e (α × V + L p ). The plastic part is defined by the motion of dislocations, both resolved and unresolved, on the current configuration and no notion of any pre-assigned reference configuration is needed.
Of significance is also the fact that MFDM involves no notion of a plastic distortion tensor and yet produces (large) permanent deformation.

Appendix B A brief review of prior work
Here, we briefly review some of the vast literature on strain gradient plasticity theories. An exhaustive review of the subject is beyond the scope of this paper.
One class of continuum models produce size-effects by modifying the hardening law to take into account the hardening due to the presence of incompatibility in elastic distortion [AB00a, BNVdG01, ATSB04, AB96, TCAS04, TAS05]. These elastic incompatibilities correspond to 'Geometrically Necessary Dislocations' (GNDs) [Ash70]. The advantage of these 'lower order gradient' theories is that the classical structure of the underlying incremental boundaryvalue-problem (bvp) with its natural boundary conditions remains unchanged, retaining all associated results on the uniqueness of solutions to the bvp of incremental equilibrium for the rate-independent (and rate-dependent) material [Hil58], as observed in [AS95]. This is in contrast with many strain gradient theories which involve higher order stresses or additional boundary conditions. The second class of models [FMAH94,FH01,NH03,Gud04] incorporate hardening effects due to plastic strain gradients by defining measures of the effective plastic strain which depend on invariants of gradients of plastic strain. These extensions of the conventional theories involve higher order stresses and require additional boundary conditions on appropriate variables.
The third class of models are based on the framework of Gurtin [Gur00,Gur02,Gur08,GR14]. In these models, the classical force balance is supplemented with a statement of microscopic force balance, one for each slip system, which also form the nonlocal flow rules of the theory when combined with thermodynamically consistent constitutive equations.
The fourth class of models [APBB04, GBB06, YGVdG04, KT08, MRR06, PHG19] augment conventional elasto-plastic theories with additional equations that aim to model dislocation motion and evolution in the domain, resulting in plastic flow. This is in contrast with the other three class of models wherein plastic flow results through constitutive assump-tions without explicitly characterizing dislocation motion and evolution. All these models [APBB04, GBB06, KT08, Gur08, MRR06, PHG19] involve the prescription of boundary conditions on scalar dislocation densities at grain boundaries in multi/poly-crystalline aggregates; this is an ambiguous, but necessary requirement, as slip system dislocation densities in different grains have different physical meanings due to the change in lattice orientation between the grains.
We now review some of the continuum models of strain gradient plasticity theory mentioned above.
Arsenlis et al. [APBB04] developed a micromechanical model of single crystal plasticity that incorporates a material length scale dependence in its framework. The authors developed a set of pdes for the evolution of slip system dislocation densities, one for each system. This requires constitutive assumptions for dislocation velocity, cut-off radii for annihilation, and segment-length interaction matrices for each component and slip system. Evolution equations including SSDs are also proposed. The framework is applied to predict size effects in simple shear for thin films using constrained and quasi-free boundary conditions for an idealized crystalline geometry.
Geers et al. [GBB06] present a strain-gradient crystal plasticity theory that has additional equations which govern the evolution of GNDs and SSDs in the domain. The back stresses due to the presence of dislocations are accounted for by using analytical integral expressions for stress fields in an infinite medium. Aspects of modeling backstresses by such assumptions are assessed in [Ach08].
Kuroda et al. [KT08] present a finite deformation higher order strain gradient crystal plasticity formulation without introducing higher order stresses. The conventional theory is supplemented by a length scale dependent back-stress which in-turn depends on the gradients of scalar edge and screw GND densities whose evolution is governed by their respective pdes. Po et al. [PHG19] present a finite deformation continuum dislocation-based plasticity formulation to model dislocation microstructure that is observed in the wedge micro-indentation experiment of [KSO + 10].
A finite strain generalization of the strain gradient plasticity formulation of Fleck et al. [FH01] is proposed and implemented in a finite element framework by Niordson et al. [NR04,NT05]. Numerically, the formulation requires the equivalent plastic strain rate to be one of the nodal variables which requires additional boundary conditions to be applied at the boundary separating parts of the body loading elastically and plastically, along with boundary conditions in higher order stresses. The model is used to study necking under plane strain tension and compression [NT05], and size effects in plane strain necking of thin sheets [NR04].
[YVdGG04] propose a nonlocal version of crystal plasticity motivated by elementary statistical-mechanics descriptions of collective behavior of dislocations [GB99,Gro97] wherein 2 coupled pdes govern the evolution of the scalar dislocation density and GNDs. The model is applied to study the shearing of a model composite material in single slip with constrained boundary conditions [YVdGG04] and bending of a single-crystal strip in plane strain [YVdGG04]. The model and computational implementation is limited to single slip deformations in 2-d problems.
Gurtin [Gur08] develops a finite deformation gradient theory of crystal plasticity involving configurational microforce balances for each slipsystem [Gur02]. A central ingredient of the theory involves inclusion of an additional term based on GND density in the free energy that characterizes an energetic hardening mechanism associated with the accumulation of GNDs.
Rudraraju et al. [RVdVG14] present stress fields of dislocations using Toupin's gradient elasticity theory [Tou62] at finite strains wherein the defects are represented by force dipole distributions. The theory involves 4 th order derivatives of displacement and therefore requires the use of isogeometric analysis [HCB05,CHB09] for its implementation. The representation of a complex network of line defects can be onerous by dipole forces. Moreover, it is not very straightforward to postulate evolution of defects coupled to the underlying stress field, given its representation in the form of dipole forces.
A hybrid method to model dislocation core widths and their mutual interactions under quasi-static and dynamic loadings was developed in [Den04,Den07]. The approach involves a combination of Pierls-Nabarro model and Galerkin methods. However, finite deformation frameworks for these approaches have not been demonstrated.