Revisiting the Application of Field Dislocation and Disclination Mechanics to Grain Boundaries

: We review the mechanical theory of dislocation and disclination density ﬁelds and its application to grain boundary modeling. The theory accounts for the incompatibility of the elastic strain and curvature tensors due to the presence of dislocations and disclinations. The free energy density is assumed to be quadratic in elastic strain and curvature and has nonlocal character. The balance of loads in the body is described by higher-order equations using the work-conjugates of the strain and curvature tensors, i.e., the stress and couple-stress tensors. Conservation statements for the translational and rotational discontinuities provide a dynamic framework for dislocation and disclination motion in terms of transport relationships. Plasticity of the body is therefore viewed as being mediated by both dislocation and disclination motion. The driving forces for these motions are identiﬁed from the mechanical dissipation, which provides guidelines for the admissible constitutive relations. On this basis, the theory is expressed as a set of partial differential equations where the unknowns are the material displacement and the dislocation and disclination density ﬁelds. The theory is applied in cases where rotational defects matter in the structure and deformation of the body, such as grain boundaries in polycrystals and grain boundary-mediated plasticity. Characteristic examples are provided for the grain boundary structure in terms of periodic arrays of disclination dipoles and for grain boundary migration under applied shear.


Introduction
Disclinations and dislocations were simultaneously introduced by Volterra, more than a century ago [1]. Dislocations are crystal defects arising from translational lattice incompatibility, as measured by the Burgers vector. Disclinations are defects originating in the rotational incompatibility of the lattice characterized by the Frank vector [2]. Disclinations have long been neglected in the field theory of crystal defects, due to their diverging elastic energy, which precludes their occurrence as single defects [3]. However self-screened configurations, such as disclination dipoles, involve bounded elastic energy and may appear frequently [4,5]. Indeed, they have been proposed for describing grain boundaries in various polycrystals, by deriving the disclination density field from crystal orientation maps, as obtained by electron backscattered diffraction [6,7]. Disclination dipoles were also exhibited in planar configurations formed by C 60 molecular layers on WO 2 /W(110) substrates [8]. Generally, disclinations are liable for complementing dislocations in the description of the lattice structure when single-valued elastic rotation fields do not exist, as in polycrystals. Disclination dipoles were shown to provide a good description of grain boundaries [9][10][11][12].
In addition to obstacles to plastic deformation by dislocation glide and twinning, grain boundaries can also contribute to plastic deformation. In ultra-fine-grained metals or in materials lacking dislocation slip systems, they can become plasticity mediators, not only by supplying dislocation sources, but also via migration mechanisms. For example, when a shear stress is applied to a bicrystal, one crystal can grow at the expense of the other by the normal motion of the boundary, inducing plastic shear in the crystal traversed by the boundary. Such shear/migration coupling is characterized by the ratio β of the displacements parallel and normal to the boundary. Surface-dislocation-based models [13,14] captured well the variations of this coupling factor as a function of the misorientation from the Burgers vector's content of the interface. Molecular dynamics methods [14][15][16] provided elementary mechanisms for boundary migration under stress. As considered here, field models involving dislocation and disclination densities are also candidates for the modeling of grain boundaries and their plasticity. The size of the investigated microstructures can be of the order of (but can be larger than) the dimensions typically considered in atomistic simulations. The field models allow analyzing the fine structure of the grain boundaries and the distribution of elastic energy at an atomic resolution length scale. It is also possible to model boundary migration by the motion of arrays of disclination dipoles, and to model dissipative and diffusive phenomena using disclination and dislocation transport equations [17].
In the present paper, we provide a review of the mechanical theory of dislocation and disclination fields [18] and show its interest in dealing with crystalline bodies subject to grain boundary-mediated plasticity. We limit the presentation to infinitesimal transformations. A more complete presentation of the model may be found in the framework of generalized disclinations [19]. By a field theory, we mean a theory of continuously distributed crystal defect densities using the tools of the mathematical theory of partial differential equations and boundary problem solving. As such, we capitalize on the earlier elasto-static theories of dislocation fields [20,21] and of dislocation and disclination fields [2]. Both theories are strongly connected: the latter reduces to the former when the disclination density field vanishes. We also benefit from the mechanical theory of dislocation fields [22], which additionally deals with dislocation-mediated plasticity. In the latter, the material displacement and dislocation density fields are derived from a set of partial differential equations complemented with initial and boundary conditions, provided constitutive information for the elastic behavior and dislocation motion is supplied [23]. This theory was shown to capture well complex features of elasto-plastic deformation such as size effects [24,25], directional hardening [24,26], plastic strain localization [27], dislocation cores [28,29] and patterning of dislocation ensembles [30,31].
When the dislocation velocity vanishes, the mechanical theory of dislocation fields [22] reduces to the elasto-static theory of dislocations [20,21]. The present framework of disclination and dislocation fields reduces to the mechanics of dislocation fields when the disclination density is ignored. It also reduces to the elasto-static crystal theory of dislocations and disclinations [2] when the velocities of dislocations and disclinations are ignored. A previous theory of plasticity accounting for both the translational and rotational aspects of lattice incompatibility was proposed by Kossecka and deWit [32,33]. Driving forces for disclination motion were discussed in [34]. The theory offered kinematic guidelines, but it was not used to solve realistic boundary value problems. In the present work, we try to use the field disclination and dislocation mechanics model in situations where the dynamic interplay between dislocations and disclinations is essential to the understanding of the plastic deformation of the material.
The outline of the paper is the following. Mathematical notations are provided in Section 2. A review of the incompatible elastic dislocation and disclination model [2] is given in Section 3. Static coupling between dislocations and disclinations occurs at this level, due to the continuity condition for dislocation densities. The equations of balance of momentum and moment of momentum are given in Section 4, which further includes the presentation of nonlinear nonlocal elastic constitutive laws that can be of interest for core problems. Section 5 is devoted to the transport properties of dislocations and disclinations. Dynamic coupling occurs at this stage, because disclination mobility yields a source/sink term for dislocations. In Section 6, guidance of the Clausius-Duhem inequality allows defining the driving forces for dislocation and disclination motion, and helps providing appropriate constitutive relations between the driving forces and defects velocities. Section 7 shows possible algorithms used for the solution of the model equations. A plane "edge-wedge" model is detailed in Section 8, and used in Sections 9 and 10 to investigate the structure and mobility of grain boundaries seen as periodic arrays of disclination dipoles.

Mathematical Notations
A bold symbol is used for a tensor, as in: A. To avoid possible ambiguities, an upper arrow can be used for a vector: V. The transpose of a tensor A is A t . Tensor subscript indices are expressed with respect to a basis (e i , i = 1, 2, 3) of a rectangular Cartesian coordinate system. Arrays of one/two dots denote contraction of the respective number of "adjacent" indices on two immediately neighboring tensors. For instance, the tensor A.B with components A ik B kj comes from the dot product of tensors A and B, and A : B = A ij B ij is the inner product. The cross product of a second-order tensor A and a vector V, the div and curl operators acting on second-order tensors are defined by: where e jkl = e j · (e k × e l ) is a component of the Levi-Civita tensor X. It is 1 if the jkl permutation is even, −1 if it is odd, 0 otherwise. The comma followed by a component index denotes a spatial derivative with respect to the Cartesian coordinate. A vector A is associated with tensor A through the inner product of A with the tensor X: The trace of the second-order tensor A is tr(A) = A ii . The symmetric/skew-symmetric parts are A sym and A skew . The hydrostatic part is A hyd = 1 3 tr(A)I, where I is the identity tensor, and its deviatoric part: A dev = A − A hyd . An upper dot represents a material time derivative.

Incompatibility in the Dislocation Model
In the proposed framework, continuity of the displacement vector field u, along with continuity of its derivatives, is assumed at any point, including between atoms, in a simply connected body that is elasto-plastically deformed, except perhaps along interfaces where continuity of the derivatives may not hold. The rotation vector is also defined at any point. These statements assume in a dual way that the material is capable of transmitting stresses and couple stresses, possibly below inter-atomic length scales. Such continuity was very early proposed [35,36] when defining pressure and stresses in quantum systems. It was later used by [37] when defining the stress tensor in quantum mechanics. By definition, couple stresses (moments of stress) at a point come from mechanical stresses existing at distant material points. As such, couple stresses introduce nonlocality in the mechanical description of the material. The distortion tensor is the gradient of the displacement U = grad u. It is curl-free: If dislocations are present, the elastic U e and plastic U p components of U contain incompatible parts, which are not curl-free, but are div-free. Indeed, if dislocations thread a patch S in the body, a constant discontinuity b in the elastic displacement exists across S, which manifests itself as a closure defect along the circuit C surrounding S: (8) b is referred to as the resulting Burgers vector of this dislocation ensemble. Applying Stoke's theorem yields: where n is the unit normal to S. The discontinuity b is also characterized in a pointwise continuous manner by introducing Nye's dislocation density tensor α: Thus, dislocation lines threading S and having a non-zero resulting Burgers vector induce the existence of an incompatible part U ⊥ e of the elastic distortion tensor. Applying the Stokes-Helmholtz decomposition of square-integrable tensor fields with square-integrable first order derivatives [38], we can uniquely define the square-integrable tensor and vector χ and z such that the elastic distortion U e reads: U e = curl χ + grad z.
The curl of U e extracts curl χ and removes grad z while, conversely, taking the divergence extracts grad z and removes curl χ. Thus, Equation (10) involves only curl χ, and curl U ⊥ e = curl curl χ = α. (13) The term grad z is the compatible (curl-free) part U e of the elastic distortion U e , and z is the compatible elastic displacement u e , up to a constant. However, Equation (13) is insufficient to identify U ⊥ e from a given dislocation field α. To ensure completeness of the identification, U ⊥ e must vanish if α = 0. Following [38,39], we thus augment Equation (13) with the supplementary condition div U ⊥ e = 0 together with U ⊥ e · n = 0 on the external boundary with unit normal n. This allows ensuring that the solution does not contain any gradient part. Applying the curl operator to Equation (13), we have: Then, using the side condition div U ⊥ e = 0, we get: Equation (15) is a Poisson-type equation for the unique determination of the unknown tensor U ⊥ e . We could have demonstrated in a similar way the existence of an incompatible part of the plastic distortion, U ⊥ p , which is opposite to the incompatible elastic distortion U ⊥ e to ensure the compatibility of the distortion U and the continuity of the displacement u. An additional compatible component U p may also exist, together with a compatible plastic displacement u p . We can then write the following equations: From the above, we see that: similarly to Equation (13). The continuity condition: derives naturally from Equations (10) and (13). The divergence-free character of α reflects the fact that dislocation lines cannot end within the material. They should either exit the material, or be closed lines. As will be shown later in Section 3.2, dislocation lines may possibly end on disclination lines. Now, in view of defining disclinations in the next section, we introduce curvatures and examine how they appear in the present dislocation model. The strain tensor is the symmetric part of U, the rotation tensor ω its skew-symmetric part, corresponding to the rotation vector ω through the relation: With this, Equation (7) yields: where grad ω is the curvature tensor κ: The trace of Equation (23) yields: div( ω) = 0, (25) and implies that tr(κ) = 0, meaning that the curvature tensor is deviatoric. Please note that transposing Equation (23) and then taking the curl yields: This is the Saint-Venant compatibility condition. Applying the above procedure to Equation (10), we find: curl e + div( ω e )I − grad t ω e = α (27) and: The incompatibility tensor η is a possible measure of the incompatibility due to the presence of dislocations. To terminate this section, we introduce the elastic curvature tensor κ e as the gradient of the elastic rotation ω e : which reflects both the continuity of the rotation vector and the curl-free character of the tensor κ e : curl κ e = 0.
Using κ e , Equation (27) can be expressed as: It is important to note for what follows that the above equation holds even if κ e is not a gradient. Finally, relations similar to Equations (29) and (31) can be obtained by using Equation (19) for the plastic strain p , rotation ω p and curvature κ p :

Incompatibility in the Dislocation and Disclination Model
We now assume that elastic rotation discontinuities ω e exist in addition to the displacement discontinuities u e = b reflected by the dislocation densities. Thus, the elastic curvature tensor κ e is not a gradient anymore. This is the situation met along grain boundaries in polycrystals. Integrating κ e along closed circuits in the manner of Equation (9) leads to the angular closure defect φ known as the Frank vector: We introduce a tensor field θ such that: curl κ e = θ (35) and φ = S θ · ndS. (36) θ is the disclination density tensor [2]. In the presence of disclinations, the elastic displacement discontinuity reads: Thus, the discontinuity u e derives from the motion of a rigid body with translation b and rotation φ, in agreement with Weingarten's theorem [40]. The Burgers vector can also be expressed as a function of dislocation and disclination density tensors: When both dislocations and disclinations exist, the incompatibility tensor η is: As similarly detailed above for U ⊥ e , we must ensure that the incompatible part κ ⊥ e of the curvature tensor is null if θ = 0. In this aim, Equation (35) must be replaced with: together with the side condition div κ ⊥ e = 0 and boundary condition κ ⊥ e · n = 0 on external boundaries. These conditions provide a unique solution κ ⊥ e of the Poisson equation: A compatible part κ e can also exist, but it is not involved in Equations (40) and (41). Therefore, one has the following equations: The incompatible plastic curvature is opposed to the incompatible elastic curvature to ensure that the total curvature is compatible and the material continuous. Then: The continuity equation for disclinations: derives from Equations (35) and (40). It mathematically expresses that disclination lines should not end within the body and form closed loops. The continuity equation for dislocations Equation (21) becomes: It means that dislocations may now terminate on twist-disclination lines. Please note that if the disclination density is null, Equations (31), (35), (38) and (49) reduce to Equations (11), (21), (27) and (30), and the elasto-static model of dislocations and disclinations reduces to the above dislocation model.

Field Equations
We consider describing the mechanical equilibrium of a simply connected body of volume V and external surface S containing dislocations and disclinations, and submitted to external loading. Since the disclination densities are associated with the incompatibility of the elastic curvatures, we are interested in a framework involving, in addition to stresses, the work-conjugates of the curvatures, i.e., the couple stresses. We therefore follow the footsteps of [41] to derive the mechanical balance equations. Let n be the outward unit normal to S, t be the traction vector (force per unit area) and m the moment vector (couple per unit area) acting on S. As they are non-essential for the present purposes, body forces and couples are neglected for the sake of simplicity, as well as inertia forces. The conservation of momentum and moment of momentum of the body therefore imposes the following conditions: Assuming that the tractions and moments at boundaries can be written as t = T · n and m = M · n, where T and M are respectively the stress and couple-stress tensors, we find by substituting in Equations (50) and (51) and applying the divergence theorem: Assuming sufficient continuity, the pointwise local equations are: Equations (54) and (55) are balance of momentum and moment of momentum equations where the stress tensor T is not necessarily symmetric. These two equations can be combined to form a higher-order balance equation coupling the symmetrical part of the stress tensor T sym and the deviatoric part of the couple-stress tensor M dev [41]: but involving neither the skew-symmetric part T skew of the stress tensor nor the hydrostatic part M hyd of the couple-stress tensor. T skew may be obtained from Equation (55), but M hyd is altogether absent from the mechanical balance equations and is therefore indeterminate in such a couple-stress theory [41][42][43].
A more general presentation of the theory involving not only curvatures and couple stresses but also strain gradients and their work-conjugates ("hyper-stresses" or "double-stresses") allows avoiding such indeterminacy. It features crystal defects more general than disclinations, and referred to as generalized disclinations. The differences between a couple-stress disclination theory involving only curvatures and a hyper-stress generalized disclination theory involving both curvatures and strain gradients are discussed in [44]. The mechanical dissipation in the material is the difference between the power of the applied forces and the rate of change of the stored energy: Here, v is the material velocity and ψ the specific free energy density. Using the balance Equations (54) and (55) together with the divergence theorem, we can write: Then we use ω i = −e imn ω mn /2 to substitute the rotation vector for the rotation tensor through Equation (22), and progressively derive: The Clausius-Duhem inequality imposes that the dissipation be positive, or zero for reversible processes: D ≥ 0. In Equation (63) the inner product between symmetric and skew-symmetric tensors is null. Hence the strain rate tensor˙ extracts the symmetric part of the stress tensor T, and the deviatoric curvature rate tensorκ extracts the deviatoric part of the couple-stress tensor. A free energy density function ψ is taken as: with contributions from the elastic strain e and curvature κ e . A more detailed description of this free energy will be given below in Section 4.3. For the time being, it is sufficient to note that, differentiating Equation (64) on the one hand, and considering elastic reversible processes such that D = 0 in Equation (63) on the other hand, leads to identifying the symmetric part of the stress tensor T sym and deviatoric part of the couple-stress tensor M dev with the partial derivatives of the free energy density with respect to the elastic strain and curvature: Equations (67) and (68) will be the constitutive relationships describing the elastic behavior of the material.

Solving the Incompatible Elasticity Problem
Consider again a simply connected body of volume V containing dislocations and disclinations, submitted to tractions and moments on its external surface S with outward unit normal n. Tractions are imposed on a part S t of S and rotations/displacements on a part S u : We are interested in describing the state of stress and couple-stress in the body in the presence of dislocation and disclination density fields (α, θ). The incompatible elastic strain and curvature ( ⊥ e , κ ⊥ e ) are assumed to be known from the solution of the Poisson Equations (15) and (41). They usually do not allow satisfying the balance Equation (56). Recalling the decomposition of the elastic strain and curvature in the sum of incompatible, ( ⊥ e , κ ⊥ e ), and compatible, ( e , κ e ), components, Equation (56) can be written as: f ⊥ is to be seen as a volumetric density of body forces associated with dislocations and disclinations. In Equations (72) and (73), the symmetric part of the stress tensor and the deviatoric part of the couple-stress tensor are obtained from the constitutive relations Equations (67) and (68). The solution of Equation (72) is the compatible elastic displacement field.

Elasticity in Dislocation/Disclination Cores
Elastic constitutive laws usually relate the elastic strain tensor to the stress tensor, both taken at the same material point. As such, they are referred to as "local" constitutive laws. They are valid in a statistical sense if the characteristic size of the volume element actually envisioned at the material point is large enough to include a large number of atoms in a perfectly ordered crystal. Local laws may fail at describing correctly elasticity if this characteristic size reduces to the order of a few atomic lattice spacings, particularly if the volume element involves crystal defects such as dislocations and disclinations. Since the mechanical behavior of a particular atom depends on its interactions with all neighboring atoms, whose location is strongly disturbed by the presence of these defects, defining the stress tensor in such a local manner may then become impossible. Nonlocal elastic laws are then needed. Nonlocal elastic laws are relationships using convolution integrals, where the stresses at a given material point depend on the deformation state in some neighborhood, whose extent is controlled by a nonlinear kernel function [45][46][47][48][49]. Such nonlocal elastic laws have been successful in retrieving phonon dispersion curves and regularizing elastic fields in the core region of crystal defects. Convolution integrals can be coupled with higher-order elastic laws to derive nonlocal elastic constitutive laws that are of interest in the core regions of dislocations and disclinations, while local laws may be sufficient in the non-defected regions. The occurrence of different elastic behavior in non-defected and defected regions was already assumed in grain and phase boundaries [50][51][52][53]. Following [54] and relying on Eringen-type convolution integrals applied to strains and curvatures, we propose below nonlocal laws describing the elastic behavior in the core region of dislocations and disclinations, and reducing to local laws in the perfectly ordered parts of the crystal.

Nonlocal Convolution Integrals
We introduce a free energy density functional at a point r in a body V in the form: .
The first two pointwise terms on the r.h.s. include symmetric positive-definite quadratic forms with A and C given at point r. This leads to the linear elastic laws used and discussed in [2,18,54] in the framework of models for dislocations and disclinations. Nonlocal contributions are added in the form of Eringen-type convolutions, which are the last two terms of the r.h.s. [55]. In such terms, the location vectors r are points distinct from point r. The difference vectors R = (r − r) and their norms R : = (R .R ) 1/2 are also defined. The spatial range of nonlocality is controlled by the nonlocal kernel K(R ), while the elasticity tensors B and D depend on (r, r ) and the elastic curvatures/strains are taken at points r . The proposed nonlocality is motivated by the strong coupling existing in the core region of defects between displacement/strain and rotation/curvature at distant locations r and r (see Figure 1). One can express the strain induced at point r by the curvature at point r: Note the presence of the above nonlocal term in the definition of the Burgers vector in Equation (37) above. Similarly, we can express the curvature at point r due to the strain at point r: The above nonlocal relations between strains and curvatures can be used to substitute curvatures and strains at points r for strains and curvatures at points r , respectively. This allows greatly simplifying the convolutions in the free energy functional Equation (74), in such a way that nonlocality applies only to elasticity tensors while all strains and curvatures are defined at point r (the reader is referred to ref [55] for more details): By differentiating Equation (77) with respect to the elastic strain and curvature tensors, we find the stress and couple-stress tensors: Using finally E ijkl = 1 2 (D ijkl + B klij ), we can write: We end up with relationships similar to those proposed in [54], except that the present elastic tensors (B, D, E) involve nonlocal integrals. As proposed by Eringen [47,48], we assume D = C. We also impose B ijkl = D klij , which allows ensuring that the free energy functional is symmetric positive-definite and that the elastic solution is stable. Interestingly, stability of the elastic strain and curvature fields is expected at the convoluted level from the nonlocal character of the formulation, whereas instability would occur in the limit of a purely local formulation [55].

Transport of Dislocations and Disclinations
In this Section, we are interested in the kinematics of dislocations and disclinations in the context of the continuous setting set forth in Section 3 for their representation by smooth density fields. In this aim, we take advantage of the transport framework for the spatio-temporal evolution of the defects density fields that derives from the conservation of the Burgers and Frank vectors during their motion. Being unquestionable from a kinematic point of view, this framework, referred to as transport of dislocation and disclination densities, provides a natural basis for the dynamic description of plasticity through dislocation and disclination motion. Since disclinations are mainly distributed along grain boundaries in polycrystals, it provides in particular a basis for describing the contribution to plasticity of grain boundary motion.
Consider a material surface S delimited by a closed circuit C. Let f be the disclination flux used to measure the rate of inflow into S of disclination lines, carrying along with them their corresponding Frank vectors φ, through a line element dx of curve C. Let V θ be the velocity of the disclinations with respect to the material, and l the unit vector along the disclination lines. Then the disclination tensor is: In the absence of a disclination source term, the conservation of the Frank's vector content demands that the rate of change of the Frank's vector of all disclination lines threading S be equal to the total disclination flux across C: For small transformations, the pointwise statement corresponding to Equation (84) is: whereθ represents the time derivative of the disclination density tensor. The rate of inflow of Frank vectors across the surface l × dx is: Indeed, suppose first that the disclination velocity V θ is in the plane (l, dx). Then, the disclination does not enter (or exit) surface S, and it is clear from Equation (86) that the flux f is zero as it should be. Conversely, f is always non-zero when V θ has a non-zero component normal to the plane (l, dx), i.e., when the disclination enters or exits surface S, and it changes sign according to whether the disclination enters or exits S. We can then obtain from Equation (86): Consequently, the local statement of balance Equation (85) becomes: Equation (88) is a transport law for the disclination density tensor θ. It can be understood as an evolution equation for θ, provided the disclination velocity V θ is known as a function of the stress/couple-stress state from constitutive statements. Its meaning is that through the curl term, the incompatible part of the disclination flux incrementally feeds the disclination density. Comparing Equations (46) and (88) it follows, after time differentiation of Equation (46) that the cross product θ × V θ can be identified with the plastic curvature rate tensorκ p , up to a gradient: At mesoscale, the termκ * p = gradω * p can be given the significance of a statistical plastic curvature rate, meaning that a non-vanishing plastic curvature rate may be obtained even when the disclination density θ vanishes at this scale [56]. Indeed, using space-time running averages of the disclination density tensor, disclination velocity V θ and plastic curvature rate tensorκ p over a domain of mesoscopic size, allows writing the mesoscopic plastic curvature rateκ p as: where overbars indicate averaged variables. It is seen thatκ p may be non-zero when the net disclination density vanishes at mesoscale (θ = 0), in which case it becomes: Dropping the overbars for convenience, the plastic curvature rate may be generally written as: As we shall discuss below, the compatible partκ p ofκ p increments the history-dependent compatible plastic curvature produced by the motion of dislocations and disclinations. Equation (92) accounts for plasticity processes in relation to the mobility of grain boundaries that cannot be described by a pure dislocation-based model.
In the theory of dislocations, the dislocation density tensor α satisfies a transport equation similar to Equation (88), in the absence of a source term [22,57]: Here, V α denotes the dislocation velocity, and the cross product α × V α can be identified with the plastic distortion rateU p , up to a gradient. However,U p is undefined in a theory of disclinations, because its skew-symmetric partω p has discontinuities. Hence, the dislocation transport Equation (93) does not account for curvature incompatibility and has to be modified in a disclination theory. For small transformations, a version correctly accounting for curvature incompatibility is readily obtained from the time derivative of Equation (33): Here the plastic strain rate˙ p is˙ where˙ * p is a gradient tensor acquiring significance at mesoscale, as discussed above for the mesoscopic plastic curvature rate tensorκ * p . The compatible part of the plastic strain rate,˙ p , is not involved in Equation (94) but it nevertheless increments the compatible plastic strain produced by the motion of dislocations through the lattice. Note importantly in Equation (94) the presence of a dislocation source/sink term, s θ =κ t p − tr(κ p )I, due to the mobility of disclinations. Such a term offers supplementary relaxation mechanisms [4,58], and will be shown to be useful later on in the prediction of shear-coupled boundary migration. When the dislocation/disclination velocities and mesoscopic contributions to plasticity (κ * p ,˙ * p ) are known from the stress/couple-stress state through constitutive relations, Equations (88), (92), (94) and (95) can be used as evolution equations for the disclination and dislocation densities. Please note that Equations (88) and (94) are first order hyperbolic equations, which implies that their solutions have propagative character, as experimentally verified in [31]. Furthermore, hyperbolicity of these equations as a strong impact on the algorithms devoted to their numerical solution from finite elements or spectral approximations [59][60][61], as detailed below in Section 7.3.

Constitutive Relations for Dislocation and Disclination Mobility
In this Section, we account for thermodynamic requirements on positiveness of the mechanical dissipation to formulate constitutive relations for the mobility of dislocations and disclinations by following a procedure proposed in [62]. Substituting Equations (65) and (66) in Equation (63), and using the additivity of elastic and plastic strain rates and curvature rates tensors, we find the dissipation to be: To identify driving forces associated with dislocation and disclination velocities, we write down the plastic strain and curvature rate tensors (˙ p ,κ p ) in terms of dislocation and disclination fluxes as in Equations (92) and (95) and substitute in Equation (96): which reads in component form: or also: Equivalently, Equation (99) is in intrinsic form: where we wrote: F α and F θ are the driving forces for dislocation and disclination motion. Using the dyadic notation α = b ⊗ t for the dislocation density (b is the Burgers vector, t the line vector), F α can be expressed as: F α = T sym · b × t, a form similar to the Peach-Köhler force for discrete dislocations. Similarly for the disclination density tensor, θ = φ ⊗ l (φ is Frank's vector, l the line vector), and F θ = M dev,t · φ × l. This is a Peach-Köhler-type relationship for disclinations. Thus, F θ drives grain boundary motion in the present disclination theory. It is produced by the couple stresses, the stresses being not directly involved, and is normal to the disclination line. Positiveness of the mechanical dissipation rate Equation (99) can be satisfied by choosing the simple constitutive relations: where (B α , B θ , λ , λ κ ) are positive material parameters. Equations (103)-(106) have the form of Newtonian viscous drag relationships, but (B α , B θ , λ , λ κ ) are generally not simple constants. They may involve non-linearity, or/and hardening laws reflecting for example thermally assisted obstacle overcoming through adequate constitutive statements. Please note that both the glide and climb components of the Peach-Köhler force are included in Equation (103). λ describes plasticity through both dislocation glide and climb at mesoscale. Similarly, λ κ describes plasticity induced by grain boundary motion at mesoscale and driven by M dev in Equation (106). Possible expressions for λ κ were proposed in [56] by averaging Equation (104).

Solution Algorithms
In the present Section, we propose two algorithms for the numerical solution of the model. The first version solves for the incompatible plastic strains and curvatures at each time step. The second version is a reduced, rate form of the latter, numerically faster, but where incompatible fields are not determined at each step.

Complete Algorithm
The complete model involves the following equations, gathered from the above sections: with the boundary conditions: The unknowns are the displacement field u and the dislocation/disclination density fields (α, θ). They are known at the initial time. The incompatible component κ ⊥ p of κ p is solution of Equation (107), with the side conditions mentioned in Equation (40). ⊥ p can be determined from Equation (108), with the side conditions div ⊥ p = 0, ⊥ p · n = 0 on S. Being history-dependent, κ p and p are initially set to zero without loss of generality. Furthermore, using the constitutive laws Equations (109) and (110), the balance of momentum and moment of momentum problem Equation (111) can be solved for the compatible displacement, as shown in Section 4.2. Then, the plastic strain rate and curvature rate (˙ p ,κ p ) are found from Equations (112) and (113), where Equations (116)-(119) have been used to evaluate the dislocation/disclination velocities and the mesoscopic terms. Only the compatible parts (˙ p ,κ p ) are employed for updating the plastic strain and curvature tensors. Finally, the dislocation/disclination densities are updated using Equations (114) and (115) and the procedure is iterated at the next time step.

Reduced Algorithm
Using the rate form of the model equations leads to a simpler and faster incremental algorithm. The main reason and difference is that the decomposition of plastic strains and curvatures into incompatible and compatible parts is possibly accounted for at the initial time, but not later. By taking the time derivative of Equations (109)-(111) and ignoring Equations (107) and (108), the algorithm uses the following equations: with the boundary conditions:Ṫ · n =ṫ d on S t (134) All fields (u, α, θ) are given at the initial time, as well as the plastic strain and curvature ( p , κ p ) arbitrarily set to zero without inconvenience. All fields can then be updated using the above rate equations. Clearly, the solutions obtained from these two algorithms may differ substantially because their continuity and differentiability properties are different.

Numerical Implementation
The finite element implementation couples a conventional Galerkin scheme for the solution of the equilibrium problem with a mixed Galerkin-Least Squares (GLS) scheme for the numerical approximation of the dislocation/disclination transport equations. Applying a standard Galerkin method to the convective transport equations leads to high frequency, propagating, numerical instabilities arising at sharp gradients. The Least-Squares (LS) finite element method, based on minimizing the norm of the residual, has the potential for eliminating such drawbacks, because it damps high frequency oscillations. Explicit LS methods have been considered so far [60,61]. However, damping is also a drawback of LS methods, because it leads to undue solution relaxation in the long-term. Therefore, an explicit Galerkin-Least-Squares (GLS) method involving convenient linear combinations of the Galerkin and LS algorithms has been proposed [60,61].
A constant drawback of finite element methods is that three-dimensional simulations are exceedingly demanding on computational resources. Spectral methods based on Fast Fourier Transform (FFT) algorithms were shown to be efficient at solving three-dimensional problems in periodic media and therefore appear as attractive alternatives to finite elements. However, the Gibbs oscillations they inherently produce in regions of strong spatial gradients reinforce the numerical instability of the solutions to the transport equations with respect to high frequency perturbations. Spectral low-pass filters in the Fourier space were shown to be efficient at damping these perturbations, and consequently at providing accurate and stable solutions to the transport problem [59]. FFT algorithms for the solution of the equilibrium problems are also available [55,63].

A Plane Edge-Wedge Model
As an illustration of the model and a basis for the next two Sections, we consider a bi-dimensional model with a distribution of wedge disclinations first discussed in [18]. In the orthonormal frame (e 1 , e 2 , e 3 ), we take the disclination tensor to be: θ = θ 33 e 3 ⊗ e 3 . The continuity Equation (48) then implies that the disclination density θ 33 only depends on the coordinates (x 1 , x 2 ). Thus, the incompatibility Equation (46) reduces to: and the Poisson Equation (47) We ignore the mesoscopic plastic curvature rateκ * p in Equation (92) as we want to deal with nanoscale problems. Therefore, the plastic curvature rates are: The constitutive relation Equation (104) for the disclination velocities then yields: and the plastic curvature rates read:κ while the transport equation becomes: In the dislocation transport Equation (94), the source/sink term s θ generates only edge dislocations (α 13 , α 23 ). Ignoring also the mesoscopic plastic strain rate˙ * p in Equation (95), the motion of these edge dislocations produces the plastic strain rates: and the plastic rotation rate:ω The above relations Equations (148) and (150) show that climb of edge dislocations (α 13 , α 23 ) produces the extension rates (˙ If all other dislocation densities are initially absent, the model therefore involves only α 13 and α 23 densities. The continuity condition Equation (49) then implies a planar distribution of these dislocations. The relation Equation (103) provides the dislocation velocities for climb: and glide: The corresponding plastic strain rates are: while the transport equation can be reformulated as: in the constitutive Equations (167) and (168), and introducing the resultant stresses and couple stresses in the balance Equations (165) and (166), one obtains two fourth-order partial differential equations for the material displacements (u 1 , u 2 ). As already mentioned, the boundary conditions comprise the prescription of tractions and moments, or/and displacements and rotations on the surface of the body. Please note that no boundary condition is specified on the disclination and dislocation densities.

Structure and Elastic Energy of Symmetric Tilt Boundaries
Copper bi-crystals with ∑ 37(610) and ∑ 5(310) symmetrical tilt boundaries of misorientation 18.9 • and 36.9 • were simulated in [64]. The bi-crystals are free of loads and periodic boundary conditions are imposed on the vertical boundaries in Figure 2. The characteristic size of the finite element mesh is 0.5Å. The wedge disclination density fields are shown in Figure 2, where they are superimposed on the relaxed atomic structure as predicted by atomistic calculations [14]. They consist of periodic arrays of wedge-disclination dipoles located at the edges of the structural units. Dilatations/contractions regions are associated with negative/positive disclinations. The local occurrence of large elastic strains suggests that a finite strain setting could be in order. Nevertheless, the predicted elastic fields are in good agreement with atomistic simulations and experiment. The elastic energy is essentially located in the dilatation/contraction areas in the core region of defects. In agreement with atomistic simulations, the energy is contained in the first atomic rows forming the grain boundary, essentially at structural units. Comparing the elastic fields of the ∑ 37(610) boundary and of the ∑ 5(310) boundary in Figure 2 shows that energy is more localized near the interface when the structural units are close to each other, due to elastic screening effects. Furthermore, the continuous description predicts features such as the shear deformation and curvature across/along the grain boundaries. As shown in Figure 3, there is a good agreement between quantitative estimates of the excess energy of grain boundaries using our model, and measurements from atomistic simulations and experiment [65][66][67]. Satisfactory agreement is obtained for all misorientations. Notably, the energy cusps predicted for the ∑ 5(310) and ∑ 5(210) grain boundaries originate in the self-screening of the elastic fields when the separation distance between structural units becomes small.

Shear-Coupled Migration of Symmetric Tilt Boundaries
As already suggested above, shear-induced normal grain boundary migration is a fascinating example of grain boundary-mediated plasticity. It was experimentally observed at moderate temperature in Al [68,69]. Following [17], we illustrate this mechanism by submitting the copper bicrystal featuring the ∑ 5(310) symmetric <001> tilt boundary studied above to shear parallel to the boundary. In Figure 4, the wedge-disclination dipoles array is observed to move downward. This motion generates edge dislocation dipoles α 13 (line direction normal to the figure, Burgers vector along the boundary). Due to the shear stress T 12 these edge dipoles glide along the boundary, generate the plastic shear strain field p 12 seen in Figure 4c and annihilate at mid-distance, as Figure 4a,b show. Meanwhile, the downward migration of the disclination dipoles produces a plastic tilt rotation ω p 3 shown in Figure 4d. Striking similitude is observed between the heterogeneity of the plastic shear/rotation fields presented in Figure 4 and those obtained by post-processing atomistic simulations [16,70,71]. Because the boundary moves downward while positive shear is produced, the shear-coupling factor β, which corresponds to the ratio between the parallel translation of the two crystals and the normal migration of the boundary, is negative. The value β ≈ −0.94 measured in Figure 4b is in reasonable agreement with atomistic simulations and experiments [14,68]. The predicted variations of the shear-coupling factor in the full range of misorientations are presented in Figure 3e. The model retrieves the two <001> and <110> branches found in experiments [68] and obtained by molecular dynamics simulations and surface-dislocation-based models [13,14]. To further analyze the role of disclinations in plasticity, the model was used for simulating shear-coupled migration in forsterite. More specifically, we have used a periodic wedge-disclination array in the particular case of a (011)/<100> tilt boundary with misorientation 60 • . This boundary was previously modeled at the atomic scale [72]. Figure 5 presents the disclination density field used to model the grain boundary structure [73]. Wedge-disclination densities are placed on the edges of structural units identified in [72]. The disclination arrangements result in a self-screening quadrupolar configuration, such that most of the energy is contained within the quadrupoles (see Figure 5b). The predicted excess energy of 1.3 J/m 2 is in good quantitative agreement with the value 1.28 J/m 2 measured from atomistic calculations [72]. When a shear stress is applied, the disclination array moves normal to the boundary, as shown in Figure 5. As detailed above for copper bi-crystals, the migration generates dislocation dipoles whose glide parallel to the boundary generates plastic shear within the region swept by the boundary. A more comprehensive study of the disclination structure, energy and mobility of grain boundaries as a function of their misorientation is clearly needed to eventually reach an accurate description of the rheology of olivine aggregates. EBSD-based disclination field imaging, as introduced in [6,7], applied to naturally deformed olivine aggregates should also be helpful in firmly establishing such models. The wedge-disclination density θ 33 (rad·m −2 ) is shown in color. The plastic shear generated by migration is shown with gray contours.

Concluding Remarks
The elasto-plastic model of dislocations and disclinations essentially consists of a set of partial differential equations, where the unknowns are the material displacement and dislocation/disclination density fields resolved at an atomic resolution length scale, and evolving at the time scale of the elementary dissipative mechanisms [18]. It is nonlocal in space and time in the conventional variables of compatible continuum mechanics and it provides a nonlocal incompatible generalization of the latter leading to well-posed problems in dislocation and disclination dynamics. It extends the elasto-plastic model of dislocations [22] by dealing, in addition to translational discontinuities and dislocations, with rotational discontinuities and disclinations. which allows grain boundary modeling in polycrystals. It also extends the elasto-static model of dislocations and disclinations [2], by dealing with plasticity mediated by transport of dislocations/disclinations within a consistent thermomechanical framework.
It has been sometimes believed that the assumption of continuity of the field variables is no longer valid at nanoscale and below, when the discrete atomistic nature of matter becomes apparent, and that continuum mechanics consequently fails to capture the physical phenomena at the scale of a few lattice spacings. In our opinion, the present results demonstrate instead that an incompatible, nonlocal continuum description involving fields of displacement and crystal defect densities that are smooth at this scale can be adequate for the purpose of describing such phenomena, in complement to atomistic representations. Indeed, the continuous disclination and dislocation theory proposed here was satisfactorily used to model the core structure of grain boundaries as well as the resulting elastic fields and excess energies. It was also shown to predict the complex interactions between disclinations/dislocations at work in the generation of plastic rotations/strains during shear-coupled boundary migration.
A continuous description of lattice defects may be more attractive than discrete approaches for various reasons. First, smoothness is desirable from the point of view of mathematical analysis and numerical computation, and because it allows dealing with defects' core properties. When viewed at a small scale, lattice defects and the corresponding distributions of elastic strain and energy are certainly better described by smooth density fields than by singularities. Building on this point of view, recent work [16,70,71,[74][75][76] has been devoted to the construction of field representations of the lattice deformation from the atomic displacements resulting from molecular dynamics simulations.
Secondly, such descriptions do not have to resolve atomic vibrations in the manner of atomistic simulations. The kinetic energy of atomic vibrations is time-averaged over periods of microseconds and rendered as dissipation, which results in higher computational efficiency. As a practical result, finite element or spectral FFT numerical approximations may allow considering the dynamics of crystal defect ensembles over time scales in the µs or more, under realistic loading rates and stresses.
Finally, continuous approaches of crystal defects at nanoscale offer a convenient basis for seamless coarse-graining by using appropriate averaging techniques, which allows building mesoscale continuous theories of grain boundary-mediated plasticity in grain aggregates [56]. This ability could allow investigating alternative disclination-based mechanisms for plasticity, e.g., triple junction motion or emission/absorption of dislocations at grain boundaries, in media where dislocation-based plasticity is limited, either by scarcity of the slip systems, as in olivine, or because dislocation glide is restricted, as in polynanocrystals.