GD3: generalized discrete defect dynamics

A mesoscale model is introduced to study the dynamics of material defects lying at interface junctions. The proposed framework couples the dynamics of discrete dislocation and disclination lines. Disclinations are expected to be natural defects at interface junctions; their presence serving the purpose of accommodating discontinuities in rotation fields at material interface junctions. Crystallography-based rules are proposed to describe the kinematics of disclination motion. A discrete-continuous couple-stress framework, in which discrete defect lines are introduced as plastic eigenstrains and eigencurvatures, is proposed to explicitly follow the dynamics of interfacial defects. The framework is then applied to study 101̄2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\left (10\bar {1}2\right)$\end{document} twin transverse propagation and thickening in magnesium. Focusing first on the case of a twin domain, It is shown that a disclination based representation of twin domains allows for an appropriate mechanistic description of the kinematics of shear transformations. In what concerns twin thickening, the stability of defects at twin interfaces is further studied. To this end, a 3D crater lying on a twin interface is described as a dipole of disclination loops. Upon self-relaxation, it is found that out of plane motion of disclinations followed by the nucleation of twinning dislocations can be activated; thereby showing that conservative non-planar motion of disclinations can be thermodynamically favorable; mechanism that had been postulated some 50 years ago.


Introduction
In the context of multi-scale modeling, the current performance of discrete dislocation dynamics (DDD) codes allows for a linkage between unit processes associated with dislocation motion -as informed by atomistic simulations-and their collective effects on both hardening and slip (Arsenlis et al. 2007;Bulatov et al. 1998;Devincre and Kubin 1997;Fivel and Canova 1999;Ghoniem et al. 2000;Zbib et al. 1998). With this, one can validate, select and calibrate, crystal plasticity based constitutive laws. Examples of success include the quantification of: latent hardening effects resulting from both short and long-range interactions between dislocations (Bertin et al. 2014;Bulatov et al. 2006;Devincre 2013;Franciosi et al. 1980;Madec et al. 2003;Kubin et al. 2008;Queyreau et al. 2009), the simultaneous effects of voids and selfinterstitial loops induced by irradiation on strength (Sobie et al. 2017a, b), etc. These advances constitute critical steps towards the development of statistically representative, microstructure and chemistry sensitive constitutive models in the very spirit of seamless multi-scale modeling as motivated in several seminal reviews in the field. Among the several roadblocks that remain, focus is placed here on the treatment of material interfaces (Bilby 1955;Bollmann 1970;Frank 1950;Priester 2013;Sutton and Baluffi 1995), and their effects on strength, hardening and eventually ductility (Capolungo et al. 2007;Hoagland et al. 2002;Mompiou et al. 2009;Wang and Ma 2004;Bayerschen et al. 2016;Spearot and Sangid 2014). One notes that in the case of dislocations, albeit with limitations pertaining to the treatment of the system dynamics that disregards entropic effects, one can trace an explicit linkage between early static and two dimensional descriptions of dislocations (i.e. line tension, stress field, effect on diffraction line broadening), molecular statics and dynamics studies of dislocation core structures and their mobility -often aided by density functional theory and minimum energy pathway search methods -and discrete dislocation dynamics. However, such linkage remains lacking in the case of defects present at material interfaces. Indeed, early as well as more recent topological descriptions of admissible material defects provide for a two-dimensional registry of likely interfacial defects (ID) in terms of disconnection, disclination or dispiration content. Let us note that other ID based material interface representations have been proposed in the literature. Some of these account for the fine scale atomic structure of the interface while other -such as those in the lineage of the Frank Bilby and Bollman's O-lattice theoryfocus mostly on the material macroscopic interface degrees of freedom. These viewpoints have largely been enhanced by further considering the relationship between the proposed description of the interface and the resulting energy. Further, the appropriateness of these constructs can be validated against atomistic simulations. This was done for example in the case of hetero-interfaces which have been recently represented as networks of discrete dislocation lines (Vattré and Demkowicz 2015).
Despite these advances, there is a clear disconnect between these interface constructs and crystal plasticity constitutive models in which interface mediated plasticity is accounted for. To address these shortcomings and to depart from a simple treatment of material interfaces as impenetrable obstacles, different classes of DDD models with increasing level of complexity have been developed. A first class of models consists of introducing dislocation sources in both bulk and material interfaces, which allows modeling the decomposition of dislocations at grain boundaries and the associated plastic straining by glide (grain boundary sliding) and climb (diffusion), as well as slip transfer (Ahmed and Hartmaier 2011;Quek et al. 2014Quek et al. , 2016. A second class of models proposes to add slip transmission and re-emission criteria at material interfaces, which can invoke critical shear stress (Fan et al. 2015), interfacial energy ) -possibly as obtained from atomistic simulations (Zheng et al. 2017)-, or line-tension (de Koning et al. 2002Zhou and LeSar 2012). A third class of models consists of explicitly modeling ID as well as their evolution and interactions with bulk dislocations. Discrete dislocation networks were recently used to model the migration of low angle mixed grain boundaries (Lim et al. 2012). Networks of misfit dislocations where used to study Ni-based superalloys (Gao et al. 2015). Discrete intrinsic dislocation networks were also used to simulate their complex interactions with extrinsic dislocations (Fan et al. 2016;Liu et al. 2011Liu et al. , 2012Yashiro et al. 2006). Discrete dislocation networks can be built from coupled kinematic/energetic considerations Demkowicz 2013, 2015;Vattré et al. 2014b;Vattré 2017a, b), or they can be atomistically informed (Burbery et al. 2017;Wang et al. 2014).
In general, ID can adopt the form of dislocations (Vattré et al. 2014b), disclinations Reinholz et al. 2016;Rösner et al. 2011) and disconnections (Hirth and Pond 1996;Pond 2006, 2016). Disconnections, consisting of a dislocation and a step (Hirth and Pond 1996), were experimentally reported in several serrated interfaces. They were shown to play a key role in mediating the motion of interfaces. This is particularly true in the cases of twin boundaries and of interphase boundaries (Mompiou et al. 2009;Khater et al. 2012;Rajabzadeh et al. 2013;Hirth et al. 2016). Atomistic simulations of the propagation of twins show that nucleation of serrations along the twin boundary occurs concomitantly with the motion of material interface junctions located at the twin tip (Xu et al. 2013). A body of work has suggested that these interface junctions could be represented by disclination lines (Xu et al. 2013).
Focusing now on the case of disclinations. These rotational type line defects have been found to be practical for describing material interfaces (Fressengeas et al. 2011;Romanov and Kolesnikova 2009). This is the case for example of high-angle grain boundaries. Alternatively, as previously stated disclinations have also been identified as the most likely natural defects in specific regions of materials. This includes triple lines, twin tips, etc. (Beausir and Fressengeas 2013;Bozhko et al. 2014;Cordier et al. 2014;Kolesnikova et al. 2016;Murayama et al. 2002;Rösner et al. 2011;Sun et al. 2016). Overall, whether disclinations are a convenient mathematical object to describe material interfaces or are intrinsic interfacial defects (similarly to dislocations) remains debated. One notes that such question is complex as the diffraction signature of disclination fields cannot be uniquely identified and as extinction conditions for disclinations cannot be derived for transmission electron microscopy. Further, while some have used disclinations as discrete static objects, the detailed crystallography of disclination lines has not yet been postulated. Thus far, networks of discrete wedge disclination lines were used to predict twin evolution in NiMnGa alloys (Reinholz et al. 2016).However a general and 3D consideration of disclination motion was not proposed. This question includes that of the natural habit planes of disclinations as well as of their out-of-plane motion. Remarkably, using a continuous description of incompatibilies in the elastic strain and curvature fields associated with disclinations, it was shown that it can be kinematically admissible for a disclination to act as a source/sink of dislocation (deWit 1970). However both the thermodynamics and crystallography of this kinematic requirement have not been studied.
To address these questions, the present manuscript proposes a simple crystallographic rule to identify potential habit planes of disclinations. Further, a generalized defect dynamic framework is derived to predict both in plane and out of plane motion of disclinations and its relationship to dislocation emission. The original work of deWit (1970; 1973), who layed-out the foundations for the generalized treatment of line defects, is therefore extended in this study to consider the full three-dimensional features associated with disclination dynamics. We propose a Generalized Discrete Defect Dynamics GD 3 model which combines the dynamics of discrete dislocation and of disclination lines. The approach can be perceived as an extension of a recently proposed discrete-continuous Fast Fourier Transform (FFT) DDD model (Bertin et al. 2015), to account for disclinations in the context of a higher-order continuum Field Disclination and Dislocation Mechanics (FDDM) framework (Fressengeas et al. 2011;. We apply the GD 3 model to the problem of the transverse propagation of 1012 twin domains in hexagonal close packed magnesium (Partridge 1967). Among others, it is shown that a disclination based construction of 3-dimensional twin domains can render the kinematics of the shear transformation. Following this, GD 3 is used to assess the stability, and the potential role of long prismatic-basal and pyramidal-pyramidal interfaces/steps on twin thickening. While the existence of these large steps has been reported experimentally (Liu et al. 2016), to date neither their stability nor their motion has been studied. In the context of a fully 3dimensional configuration, these defects are represented as a crater on the twin interface which is formally described with a disclination dipole arrangement. As will be shown, the simulations suggest a new mechanism whereby large steps along the twin interface could decay by the emission of twinning dislocations. The results are qualitatively similar to those predicted in Xu et al. (2013) via atomistic simulations and prove that the aforementioned dislocation source mechanism, proposed by DeWit (1970) is not only kinematically admissible but also thermodynamically possible.

Generalized discrete defect dynamics
Useful mathematical notations are provided in the appendix. A small strain setting is assumed. The material displacement field vector u is single-valued and continuous at any point of the crystal. The distortion tensor, defined as the gradient of the displacement U = grad u, is curl-free: (1) The strain tensor is the symmetric part of U, the rotation tensor ω its skew-symmetric part and the associated rotation vector ω is: It is also a single-valued continuous field. The gradient of the rotation vector is the second order curvature tensor κ: The total curvature tensor is a compatible curl-free tensor, like the distortion tensor. The total distortion and curvature tensors are additively decomposed into elastic and plastic parts: By decomposing the elastic and plastic distortion into strains and rotations, one also has the relations: In the above two equations, the right part corresponds to the the curl of the elastic and plastic spin tensors, respectively.
The mechanical framework proposed considers line defects (i.e. dislocations and disclinations) both as discrete and as continuous objects. As such the mechanical fields (i.e. Cauchy stress, couple-stress, total strains, total curvatures) are obtained by first quantifying the plastic strain and plastic curvatures associated with the motion of dislocations and disclinations, respectively. In the following, when needed, the symbols and ⊥ will be used to differentiate between disclination and dislocation related quantities, respectively. In the case of dislocations lying on specific slips systems, one writes the total plastic distortion field within the medium as follows: The symmetric part of the plastic distortion is the plastic strain while the skew-symmetric part is the plastic spin In the above relations, b s and n s⊥ denote the unit Burgers vector and slip normal direction of slip system s. The local plastic shear γ s is obtained by multiplying the Burgers vector of the dislocation considered by the area is has sheared, normalized by the simulation volume. Importantly, Eq. (10) describes the crystal lattice rotation due to dislocation glide. Lattice rotations are therefore accounted for as they will impact the local elastic fields as a consequence of the elastic anisotropy of the crystal. Further, the evolution of lattice curvatures (the gradient of lattice rotations) due to dislocation glide will also affect the dynamics of disclinations.
We decompose the plastic curvature into disclination and dislocation parts: The plastic curvature tensor due to dislocations is fully compatible as it is the gradient of the plastic rotation. It can be written as: On the contrary, the plastic curvature due to disclinations κ p must contain both an incompatible part κ p i and a compatible part κ p c . The compatible part will be defined later on.
The net plastic rotation jump due to the presence of a disclination is denoted with the Frank vector . Assuming known the planes on which disclinations can move, similarly to the case of dislocations, one can thus write the local incompatible plastic curvature field as follows: where s and n s denote the unit Frank vector and 'slip' normal direction of system s. Obviously the 'slip' planes for disclinations are not necessarily the same as those of dislocations. The plastic curvature κ s is obtained by multiplying the Frank vector of the disclination by the area it has swept, normalized by the simulation volume. In doing so, we impose that disclination lines and loops move in crystallographic planes, which is physically realistic and original as compared to the Field Disclination and Dislocation Mechanics continuous framework (Fressengeas et al. 2011). As shown in the next section, it is important to note that the plastic curvature induced by disclinations also generates plastic strains, which in turn affect the Cauchy stress field. Given plastic strain and curvature fields due to a distribution of dislocation and disclination lines, we now search for the associated internal stress and couple-stress fields. In the absence of significant body forces, both momentum and moment of momentum are conserved (Cosserat and Cosserat 1909): where the stress tensor T is generally non-symmetric. Thus, the rotated stress vector T : X balances the divergence of the couple-stress tensor M when T is non-symmetric. Given a description of the specific free energy density function containing contributions from elastic strains and curvatures, one can introduce the constitutive relationship for the symmetric Cauchy stress T sym and deviatoric couple-stress M dev ): Plastic strains and curvatures due to defects appear in the above balance equations, through the additive decomposition e = ( − p ) and κ e = (κ − κ p ). Further note that tensors A, B and D involve internal characteristic lengths that are the order of a fraction of the Burgers vector magnitude Po et al. 2014;Seif et al. 2015).
Focus is now placed on the motion of discrete dislocation and disclination lines under the action of internal stress and couple-stress fields. At any point on the dislocation, the Peach-Koehler force is written as: where b is the Burgers vector and t ⊥ denotes the line unit vector. Similarly, the Peach-Koehler-type force acting on a point lying on a disclination is expressed as follows: where is the Frank vector and t denotes the line unit vector. It is a discrete version of the driving force defined in earlier work (Fressengeas et al. 2011). Naturally, in the context of discrete defect dynamics, forces (18) and (19) should be integrated along a discretized segment in order to extract the forces acting on a node. From the known Peach-Koehler forces, the local velocity of points lying on the dislocation or disclination lines can be computed. In the present case, one assumes an overdamped equation of motion expressed as follows: where B ⊥ and B respectively denote the viscous drag coefficient matrix for dislocation and disclination motion, while V ⊥ and V are the defect velocity vectors. Atomistic simulations can be used to estimate the drag coefficients for dislocation glide, depending on temperature, strain rate, segment character etc. (Bitzek and Gumbsch 2004;Cho et al. 2017). However, they are not known for disclinations. We will assume in our simulations that dislocation mobility is significantly higher than disclination mobility, meaning that disclination-mediated interface motion and plastic relaxation mechanisms is then much slower than dislocation glide in grains. Further for the sake of simplicity, the mobility of dislocations and of disclination are taken as independent of the segment character. Here as well, velocities (20) and (21) are then rewritten by introducing the interpolation functions used to discretize each segment. The resulting equations are then integrated over each segment in contact with a specific node such as to extract an equation of motion at each node as is typically done in discrete dislocation dynamics.
Finally, from the viewpoint of kinematics alone, one can introduce a generalized Burgers vector conservation law. It simply states that the flow of Burgers vector is a conserved quantity. This type of expression is typically used to solve the kinematics of junction formation and unzipping processes, i.e. b i = 0. However, in the present case the equation of continuity is expressed dynamically and in the presence of disclinations, such that non-curl-free plastic curvature rates can lead to the generation of new dislocation lines. This was originally proposed by deWit (1970) and recently derived in the Field Dislocation and Disclination Mechanics model (Fressengeas et al. 2011). In order to express this dislocation source/sink term, we first define the Nye dislocation density tensor component α ij = 1 S nb i t ⊥ j in the cartesian frame. The number n of dislocation lines of Burgers vector component b i gives a total length of Burgers vector per unit resolution surface S. With this areal tensor, the evolution of dislocation densities writes: Further, recalling the decomposition of the plastic distortion into strains and spins, one has:α The plastic curvature rates include contributions from dislocation glide (the gradient of the plastic rotation) and from disclination motion, such thatκ p =κ p⊥ +κ p . The last contribution is the source/sink term for dislocations (Fressengeas et al. 2011). Coming back to the discrete version, this source/sink term can be written in component form as: In the above discrete form,ṅ is the number of generated dislocations. The vectors b i and t j are expressed in the reference coordinate system. As such, they have to be projected on possible crystallographic dislocation systems to nucleate discrete dislocations. A priori, there is no unique solution (nucleation of dislocation partials for instance) and we shall select the solution, based on energetic arguments such as dissipation, or on experimental observations and atomistic simulations. A first application of this source/sink term is shown later on in this paper.

Numerical implementation
In this section, we present the important numerical aspects related to the mesoscale model. The field equations will be solved by using Fast Fourier Transform (FFT) algorithms. First, we present the regularization method used to distribute the plastic fields induced by discrete dislocation and disclination lines (Bertin et al. 2015) within the FFT grid. Second, we show how to uniquely calculate the plastic strains arising from plastic curvatures induced by disclinations (Fressengeas et al. 2011). Third, we provide details on the spectral FFT code used to solve the balance equations.
We start with the calculation of plastic eigendistortions U p⊥ and plastic eigencurvatures κ p due to discrete dislocation and disclination lines, respectively. The reader is referred to the work by Bertin et al for more details (Bertin et al. 2015). The numerical method is presented in the case of dislocation lines, then the exact same procedure applies for disclination lines. Nabarro showed that dislocations can be described as a set of coherent misfitting platelet inclusions (see for instance Mura (1987); Wang et al. (2001)). This representation of dislocations was recently used to develop the so-called Discrete Continuous Model (DCM) (Lemarchand et al. 2001;Vattré et al. 2014a). In this eigenstrain-based numerical formulation, the elastic fields of dislocations is obtained by introducing dislocations as plate-like inclusions (see Fig. 1), bearing the appropriate plastic strain, in a finite element solver. With this approach, a dislocation line -defined by its Burgers vector b and its defect surface S with unit normal n -produces a plastic distortion: where δ(S − x) denotes the three-dimensional Dirac delta function that is zero everywhere except on surface S, and which accounts for the displacement discontinuity [u] = b across S. The plastic distortion U p (x), directly resulting from the motion of dislocation lines, is numerically updated using the same regularization procedure. The increment of plastic distortion generated at voxel x d from the motion of dislocation segments on all slip systems is expressed as: where the summation is performed over all slip systems s with Burgers vector b s and unit normal n s , and where dγ s (x d ) denotes the plastic shear increment resulting from crystallographic slip on system s. Considering the motion of all dislocation segments across voxel x d , increment dγ s (x d ) is expressed as:  (27) where dγ x d ij denotes the shear produced by the glide of dislocation segment ij at voxel x d , and where the summation is performed over all segments ij gliding across voxel x d . As depicted in Fig. 1, dislocations are represented as a plate-like inclusion of thickness h = 2L, with L denoting the distance between voxels. With this, the shear produced by the motion of a dislocation can be determined as follows. It is considered that an elementary sheared area dS(x) centered in x swept by the glide of a portion of a dislocation segment ij produces an elementary homogeneous plastic shear dγ (x) within an elementary spherical volume dφ(x) of radius h/2 (Vattré et al. 2014a), such that: where b denotes the magnitude of the Burgers vector of the dislocation line. c(d) denotes a correction factor, as defined in Bertin et al. (2015), which ensures that the predicted stress fields remain correct regardless of the distance d between the segment and the center of the voxel. Following Eq. 28, the plastic shear dγ x d ij in expression (27) produced by the glide of a dislocation segment ij is regularized at each grid point x d as: where V e = πh 3 /6 is the volume of the elementary spherical sheared volume dφ(x d ) of radius h/2, and where quantity dS x d ij corresponds to the intersection between the area dA ij swept by dislocation segment ij and the elementary sphere of volume V e centered at voxel x d . Further details on this regularization procedure (e.g. choice for the regularization parameter h, numerical procedure to compute dS x d ij , numerical implementation, etc.) can be found in a recent paper (Bertin et al. 2015). The same procedure applies for disclinations, the Frank vector substituting for the Burgers vector and using disclination slip systems.
From the known incompatible plastic curvature due to disclinations κ p i , one must quantify the associated incompatible plastic strain p and compatible plastic curvature κ p c fields. The reader is referred to the original paper introducing Field Disclination and Dislocation Mechanics for more details (Fressengeas et al. 2011). Once p is known, the total plastic strain p = p + p⊥ is obtained and the balance equations can be finally solved. The plastic strain p satisfies: The compatible curvature κ p c is the gradient of a continuous plastic rotation field. The latter rotation and the plastic strain p compose the plastic distortion U p , such that the above incompatibility equation can be written as: Note that U p cannot contain any gradient part to ensure uniqueness of the solution (Acharya 2001). As such, it must be divergence-free. Hence, by taking the curl of above equation, U p is solution of the higher-order equation: div gradU p = U p = −curl κ p i t − tr κ p i I .
Note in the above equation that the symbol in front of the plastic distortion indicates the Laplace operator applied to the plastic distortion. No boundary condition on U p applies here because of periodic boundary conditions in our FFT framework. Equation (32) is solved in the Fourier space by using FFT algorithms and special rules for evaluating spatial derivatives. The numerical method used is presented in a recent paper (Berbenni et al. 2014). Finally, the plastic strain p can be derived as the symmetric part of U p . The compatible plastic curvature κ p c can be calculated from the skew-symmetric part of the distortion.
Once the above plastic strains and curvatures are known, we can solve for the stress and couple-stress balance Eqs. (14,15). Recall that the elastic moduli tensors A, B and D involve internal characteristic lengths that are the order of a fraction of the Burgers vector magnitude. Hence, at a resolution scale sufficiently larger than interatomic distances, which is the case in the proposed mesoscale model, the couple-stresses do not affect the balance of the stress field and the skew-symmetric stress is negligible as compared to the symmetric Cauchy stress. Hence, solving the above balance equations or simply div(C : e ) = 0 yields exactly the same stress and couple-stress fields. Obviously, it is not true anymore when using a very small resolution scale to model the core elastic fields of defects (Taupin et al. 2017). In the model, the elastic stiffness tensor C is heterogeneous because of crystal orientation and elastic anisotropy. We use the accelerated FFT spectral scheme formerly introduced by Eyre and Milton (Eyre and Milton 1999;Michel et al. 2001) and recently used in the DCM-FFT-DDD code (Bertin et al. 2015). The converged solution of the balance equations provides the total strain tensor .
The last step consists of deriving the total curvature tensor κ, such that the elastic, total minus plastic, curvature is known and the stresses and couple-stresses can be evaluated as per constitutive elastic laws (16, 17). The strategy proposed to derive the curvature tensor in the Fourier space is as follows. The Cauchy stress field satisfies the balance equation: which reads in component form We now switch from the real space to the Fourier space. Letλ denote the Fourier transform of a function λ and ξ be the Fourier frequency vector with components ξ i in a reference cartesian frame. The imaginary number is denoted by i. The balance Eq. (34) becomes The heterogeneous elastic moduli tensor is C = C 0 +δC, where C 0 is the stiffness tensor of a homogeneous reference medium and δC is the spatial fluctuation of C. Expressing the Cauchy stress as T sym ij = C ijkl u k,l − p kl allows rewritting Eq. (35) in the well-known form whereτ ij is the stress polarization tensor that includes both the fluctuation δC and the plastic strain. The above equations set an algebraic linear system of the form whereĜ is associated with the homogeneous reference tensor C 0 , andû k is the material displacement. The latter is given bŷ From the expression of the displacements, we now derive the curvature components. We first express the rotation vector as half of the curl of the displacement vector: The curvatures thus read: We finally cast the stress polarization tensor inf. Eq. (40) is transformed into : By defining the following fourth-order tensor: The total curvature finally reads: An inverse FFT can finally be used to obtain the curvature tensor in the real space. Finally, one obtains the elastic curvature tensor by subtracting the plastic part. Note that the material rotation vector is also known from Eq. (39), which allows calculating the elastic rotation by subtracting the plastic rotation.

Illustrations
We now present two illustrations of the proposed model. A simple cubic crystal structure with lattice parameter a = 0.32nm is assumed. Isotropic elasticity is used, with shear modulus μ = 17GPa and Poisson coefficient ν = 0.3. Further, we consider a periodically reproduced cubic simulation box with dimension (1800a) 3 . The critical distance for discrete defect interactions (annihilation, junctions etc.) is set to be 10a = 3.2nm and the defect remeshing distance to 20a = 6.4nm. The FFT grid is made of 128 3 voxels, implying a spatial resolution for the elastoplastic fields of ∼ 14a ∼ 4.5nm. Regarding defect velocities, the drag coefficient matrix is taken as diagonal and isotropic such that only one drag coefficient is necessary to describe dislocation glide. Here we take B ⊥ = 10 −2 Pa.s. The drag coefficient matrix for disclinations is also taken as diagonal and isotropic resulting in the need for only one drag coefficient. Its value is 100 times that of dislocations. In these illustrations, we restrict ourselves to the sole consideration of disclination dynamics. The glide planes of disclinations are limited to cube faces of the crystal cell; e.g. {001}. The Frank vector of disclinations is aligned with the cube directions and its magnitude is denoted by . In the first simulation, we propose to model an isolated circular disclination loop. Figure 2 shows the self-relaxation of a circular disclination loop of radius r=100nm, under the action of its own couple-stress field. The Frank vector is along [100] direction with magnitude arbitrarily chosen as = 1rad. The initial loop is shown on the top right and top left of the figure. The loop can glide in the (001) plane. The evolution of the shear stress T 23 (respectively the couple-stress M 12 ) during the self-relaxation of the disclination loop is also shown on the left (respectively right) column in the figure. The simulation predicts the self-annihilation of the disclination loop under the action of its own couplestress field. In short, it is shown here that the planar motion of disclinations is similar to those of dislocations. Although not shown here, applying a couple-stress M 13 leads to the expansion of the disclination loop. Focus is then placed on the case of short and long range elastic interactions between twin disclinations. Consider a two-dimensional scenario in which a matrix contains a stack of domains as depicted in Fig. 3. Further assume that the crystal disorientation axis across each domain is parallel to X 1 . Each domain is bounded by a wedge disclination quadrupole, the wedge disclination lines being infinitely straight. The Frank vector is along [100] direction with magnitude arbitrarily chosen as = 0.1rad and the elastic fields are invariant w.r.t. the X 1 axis. This configuration is very similar to that obtained in 2D simulations of twins in shape memory alloys (Reinholz et al. 2016), except that it is more symmetric here. The initial in plane internal shear stress field is shown in the Fig. 3 and resembles to that predicted in simulations of shape memory alloys. To model twin thickening/thinning, the disclination slip planes are now taken as (010). The system is then set to relax in its own elastic field. Figure 3 shows the migration of twin boundaries in the normal [001] direction, i.e. thickening or thinning of domains, until one recovers a defect free matrix. We note that the long and short range (i.e. annihilation) interactions between the twin domains are fully accounted for. Fig. 3 Self-relaxation of twin stacks. Initial twin stacks are let to self-relax by thickening/thinning of twins in the normal direction. a: sketch of the simulated twin stacks using wedge disclination lines. b: Colour-coded is shown the initial in plane shear stress field (MPa) due to disclination lines (red lines). c View in the (X 1 , X 3 ) plane showing the progressive growth and shrinkage of twins during self-relaxation until a defect-free matrix is obtained

Applications
Focus is now placed on the application of the disclination and dislocation dynamics framework to study the kinematics and thermodynamics of twin propagation in Mg, as well as twin thickening mediated by the motion of ledges and emission / propagation of twinning dislocations at twin interfaces. We consider 1012 twins (Partridge 1967). The twin plane normal is 1012 and coincides with the X 3 axis, the shear direction is 1011 and coincides with the X 2 axis. The normal to shear direction in the twin plane is 1210 and coincides with the X 1 axis. The value of twinning shear is s = 0.131. Twinning dislocations gliding in the twin plane have a Burgers vector in the shear direction 1011 and the magnitude is b = 0.0446nm. The elementary step associated with disconnections and separating adjacent twin planes is h = 0.378nm. In the reference coordinate system chosen, the shear transformation s leads to a distortion U p 23 = s. The plastic distortion can therefore be decomposed into plastic strains, p 23 = p 32 = s/2, and spins ω p 23 = −ω p 32 = s/2. In consequence the rotation vector is p 1 = −s/2. Therefore, the magnitude of the Frank vector of the elementary disclinations used to build the twins and twin interface ledges is 0.0655rad=3.7°. The magnesium lattice parameters are a = 0.32nm and c = 0.52nm. Isotropic elasticity is reasonably assumed (Capolungo et al. 2010), with shear modulus μ = 17GPa and Poisson coefficient ν = 0.3.
We postulate that disclinations could be natural defects at the intersections of material interfaces (i.e. junction of interfaces). This is the case triple lines, twin tips or ledges on the twin interface. This has also been suggested in other body of works (Bollmann 1970;Xu et al. 2013; Barrett and El Kadiri 2014). With regards to disclinations' habit planes, as suggested by atomistic simulations -albeit providing at best circumstantial evidence-we further postulate that disclinations can naturally glide on minimum energy interface planes. In the case of tensile twinning in magnesium, disclinations would therefore glide on the 1012 twin plane. In addition, wedge disclinations which would bound the tip of a twin domain could also glide/climb on the basal/prismatic interface. The latter was shown to be a minimum energy interface. In what concerns the 'dark side' of the twin (Liu et al. 2016), i.e. observing the twin domain along the twin shear direction, twist disclinations would then potentially 'cross-slip' on pyramidal/pyramidal boundaries (i.e. 2 110 // 2110 or potentially 1 010 // 1011 ).

transverse propagation of the 1012 twin
As stated in the above, disclination and their motion allow for an explicit description of lattice rotations. As such they appear as ideal defects/constructs to track the kinetics of propagation of twin domains. This is studied in the case of an isolated 1012 twin domain which stability and transverse propagation is modeled as shown in Figs. 4-6. A twin embryo is created simply by placing two disclination loops of opposite Frank vector 1 = e 1 = − p 1 = ±s/2 above one another on parallel twin planes. The positive (resp. negative) loop is placed on the bottom (resp. top) of the dipole, such as to delimit a twin domain with negative plastic rotation −s/2 with respect to the matrix. Figure 4 shows the elastic fields within the initial twin domain. In the top of the figure are shown the in plane shear stress field T 23 and the elastic dilatation/contraction in the plane parallel to the twinning shear direction 1011 . Disclination lines tangent to that plane are of pure wedge character. The location of elastic dilatations (resp. contractions) allows identifying negative (resp. positive) wedge disclinations. The bottom of Fig. 4 shows the out of plane shear stress fields in the plane normal to the shear direction. The disclination lines In blue colour is shown the plastic rotation around X 1 . Only plastic rotation values smaller than -0.06 rad are shown tangent to that plane are of pure twist character. One can see that the plastic eigenstrains and eigencurvatures resulting from the generation of the twin domain is compensated in the matrix phase such as to lead to the generation of large 'back stresses' acting on the disclinations. Importantly, one can note an asymmetry of the shear stress field T 23 in planes parallel (Fig. 4a) and normal (Fig. 4c) to the twinning shear direction. Strong positive stresses can be seen in front of the twin in the plane parallel to the shear direction, while they are much less significant in the other plane. Hence, from a purely elastostatic standpoint and considering only stresses, the latter suggest that twin propagation should be favored in the shear direction 1011 as compared to the normal direction 1210 . However, from a dynamic point of view, disclination motion is driven by couple-stresses and the kinetics of propagation of the twin domain in shear and lateral directions depend on the mobility of disclinations. As such, introducing a segment-character dependent disclination drag will allow controlling the kinetics, together with internal elastic fields of the twin. This will be the object of further studies. Figure 5 shows the initial twin domain as well as the plastic rotation field (the matrix phase is taken as reference for computing orientations). It then shows the evolution of the disclination loops and the associated rotation field during self-relaxation under the internal couple-stress field of the twin. The twin domain is consistently observed to shrink, together with the rotation field, such that one ends up with a single phase matrix. However and as expected, upon imposing a macroscopic shear strain to the system, 23 = 0.01, it is found that the twin domain can now expand. This is shown in Fig. 6. This applied shear strain yields a couple-stress as per non local elasticity tensor B (Taupin et al. , 2017, which is in turn responsible for disclination loops expansion. Very importantly, not only the prediction of twin growth in the transverse direction can be predicted by our 3D model, but so does the evolution of the rotation field, as shown in Fig. 6. We note that the shape of the twin domain remains almost perfectly circular during shrinking and growth.

Twin thickening
Thickening of a twin domain can be mediated by several processes which all lead to the generation of ID. Focus is placed on the generation of twinning dislocations from ledges present on a twin interface. For this first study, the intent is not to extract the kinetics and overall contributions of these mechanism but solely to assess the thermodynamics and kinematics of the mechanism. As shown in Fig. 7, we consider a twin plane containing a 'crater' in which the matrix phase locally penetrates into the twin domain. The boundaries of the crater can therefore be modeled as a disclination loop dipole with Frank vector corresponding to half of the twinning shear as in the previous case. For the sake of simplicity a square loop configuration is taken and the initial loop sizes are 200a and 100a. The crater height is 50a. The four interfaces connecting the two square disclination loops can be interpreted as Basal / Prismatic (BP) and Pyramidal/ Pyramidal (PP≈ PP1 or PP2) interfaces. As in the previous case, disclination loops can glide in the 1012 twin plane. However, out of plane motion of disclination loops is also allowed in the BP and PP interfaces. The choice of the plane of motion is simply based on the driving forces magnitude acting on disclination lines. If disclinations experience out of plane motion, following Eq. (24), twinning dislocation loops with Burgers vector b in the twinning shear direction will be generated in adjacent twin planes once disclination loops have moved by an elementary step distance h perpendicular to the twin plane. Upon self-relaxation, the couple-stress field induces an out of twin plane motion of the disclination loops crossing onto the BP and PP1 planes. The loops attract each other, such that the upper loop moves down while the lower loop moves up. Once both loops have moved by a distance h perpendicular to the twin plane, a dipolar arrangement of square twinning dislocation loops (red lines in the figure) is thus generated. As the mobility of dislocations is higher than that of disclinations, the upper and lower dislocation loops rapidly expand and shrink, respectively. Interestingly, the upper loop reaches an equilibrium elliptical shape, as shown in Fig. 7d, due to a balance between the repulsing stress field of the crater and the line tension of the elliptical dislocation loop. The loop is more elongated in the screw direction. As sketched in the Fig. 7, we shall interpret the observed mechanism as the nucleation of disconnection loops from ledges in the twin interface. The generated disconnection loops are indeed composed of a step and of a dislocation, and they help propagating the twin interface. Such a mechanism is very similar to observations in twin boundary interfaces in Mg as obtained from MD simulations, and can be invoked in other interfaces (Mompiou et al. 2009;Khater et al. 2012;Rajabzadeh et al. 2013;Hirth et al. 2016).

Conclusions
We have proposed a novel Generalized Discrete Defect Dynamics GD 3 mesoscale model. In addition to discrete dislocation lines, discrete disclination lines are introduced to account for the incompatibilty of elastic curvatures and discontinuity of elastic rotations in material interface defects, such as triple lines, twin tips, and interface ledges. Apart from the numerical implementation, our work proposes that disclinations could be intrinsic defects that are to be unambiguously found at junction between material interfaces. Further, our work postulates that disclination motion is constrained to be located on minimum energy interface planes (in the case of twinning the Coherent Twin boundary, Basal Prismatic boundary etc.) Third, we evidence the possibility for disclinations to change glide plane and emit dislocation in a fully conservative fashion. Applying the model to 1012 twinning in magnesium, It is shown that twin domains can be appropriately described by using dipoles of disclination loops. Further, both shrinking and propagation of twins can then be simulated. Importantly, the introduction of disclination loops as plastic eigencurvatures allows to conveniently track the evolution of rotations during twin propagation. The disclination based constructs of twins thus appears to be practical as internal stress fields, rotations and propagation can be rendered simultaneously. Further, the model was used to study the stability of 3-dimensional defects, craters, along the twin interface. Simulations suggest that large ledges along the twin boundaries could decay by generating twinning dislocations; a phenomenon observed in atomistic simulations that insofar had not been mechanistically described or understood. The prediction of this mechanism is the sole consequence of the higher order kinematics and thermodynamic framework of the model. Naturally, the proposed framework pauses several interrogations pertaining to the mobility of disclinations, their rate of decay and the short range interactions with dislocations. These will be pursued in further studies.

Appendix: mathematical notations
A bold symbol denotes a tensor. When there may be ambiguity, an arrow is superposed to represent a vector: V. The symmetric part of tensor A is denoted A sym . Its skewsymmetric and deviatoric parts are A skew and A dev respectively. The transpose of a tensor A is denoted A t . The tensor A · B, with rectangular Cartesian components A ik B kj , results from the dot product of tensors A and B, and A ⊗ B is their tensorial product, with components A ij B kl . The vector A · V, with rectangular Cartesian components A ij V j , results from the dot product of tensor A and vector V. A : represents the trace inner product of the two second order tensors A : B = A ij B ij , in rectangular Cartesian components, or the product of a higher order tensor with a second order tensor, e.g., A : B = A ijkl B kl . The cross product of a second-order tensor A and a vector V, the div and curl operations for second-order tensors are defined row by row, in analogy with the vectorial case. For any base vector e i of the reference frame: In rectangular Cartesian components: (curlA) ij = e jkl A il,k .
where e jkl is a component of the third-order alternating Levi-Civita tensor X. A vector A is associated with tensor A by using its trace inner product with tensor X: In the component representation, the spatial derivative with respect to a Cartesian coordinate is indicated by a comma followed by the component index. A superposed dot represents a material time derivative.