A new method to model dislocation self-climb dominated by core diffusion

The mobility of atoms in dislocation core regions is many orders of magnitude faster than in the surrounding lattice. This rapid atomic transport along dislocation cores plays a signiﬁcant role in the kinetics of many material processes, including low-temperature creep and post-irradiation annealing. In the present work, a ﬁnite element based analysis of the dislocation core diffusion process is presented; based on a variational principle for the evo- lution of microstructure. A dislocation self-climb model is then developed by incorporating this ﬁnite element core diffusion formulation within the nodal based three-dimensional discrete dislocation dynamics framework. The behaviour of an isolated loop in bcc iron is brieﬂy reviewed, and simulations are extended to include the loop coarsening processes of both parallel and non-parallel loops by self-climb plus glide mechanisms, in which the huge time scale separation between climb and glide is bridged by an adaptive time step- ping scheme. Excellent agreement is obtained between the numerical simulation, the theoretical solution of rigid prismatic loops and published experimental results. The coarsening process of a population of loops is simulated to investigate the mechanisms of the accumulative interactions and large-scale-patterning in bcc materials. angle for the climb of near screw dislocations, or the decou-pling of the glide and climb in each increment, the proposed methodology sheds new light on the simulation of core diffusion dominated mechanical processes and enables physically based predictions of post-irradiation annealing with high computational eﬃciency. This new method will be used in future work to investigate low-temperature creep in single crystals ( Frost and Ashby, 1982 ). Another interesting topic is the so-called displacive plasticity for micro- and nano-sized crystalline materials ( Li and Ma, 2018 ), where the strength is controlled by diffusional plasticity. We would like also to emphasise that although the present study focuses on the post-irradiation annealing process in a bcc Fe system, the underlying principle, and the self-climb/glide mechanism, are general and could be used to simulate dislocation activity and inelastic deformation in many different materials.


Introduction
Dislocations are the predominant carriers of plastic deformation in crystalline materials ( Anderson et al., 2017 ), which bridge the atomistic-scale deformation events with the macroscopic strength and plastic properties of crystalline materials. It is generally believed that dislocation glide motion, which follows the phonon-drag law ( Cai and Bulatov, 2004 ), dominates at low temperature. While, at elevated temperature, climb motion normal to the slip plane plays a critical role, where the plastic strain is produced by the glide motion at a rate controlled by climb ( Caillard and Martin, 2003 ). However, significant climb motion has been reported at intermediate temperatures which are lower than half of the melting point T m ( Mompiou et al., 2003;Somekawa et al., 2003 ), or even at room temperature ( Hay, 2004;Wang et al., 2009 ). The way in which point defects interact with the dislocation structures significantly affects the plastic behaviour of crystalline materials. It is, therefore, necessary to understand, and be able to model, these processes in an appropriate framework to predict the influence of climb on mechanical properties. Dislocation climb requires absorption or emission of point defects. The point defects are transported either along the dislocation line, in the core region, known as core (or pipe) diffusion; leading to a conservative climb motion called selfclimb ( Anderson et al., 2017 ). Or perpendicular to the dislocation line into the surrounding lattices, known as lattice (or bulk) diffusion and the corresponding motion is the more commonly studied non-conservative climb. As illustrated by  , in the dislocation core region, the strain is high and atoms are loose-packed, so that atoms in this region have a higher jump frequency ( Pun and Mishin, 2009 ) than atoms in the regular lattices, leading to a much faster diffusivity in the core region. Diffusivity in a dislocation core is reported to be six orders of magnitude greater than that in the bulk  and it might be even more rapid than can be observed by tracer techniques ( Turunen and Lindroos, 1974 ). In fact, the effective diffusivity of a solid, in which diffusion occurs by both lattice diffusion and core diffusion, can be calculated as ( Frost and Ashby, 1982 ), D eff = D l f l + D core f core (1.1) where D l is the lattice diffusion coefficient, f l is the fraction of atom sites associated with lattice diffusion, D core and f core are the corresponding quantities for the core. f core is determined by the cross-sectional area of the dislocation core a core and dislocation density ρ, that is, At steady-state, ρ is a function of stress and temperature, given as ρ = A ( σ s / Gb ) 2 by Argon (1970) , with σ s denoting the effective shear stress, G denoting the shear modulus and b is the magnitude of the Burgers vector, A ≈ 1, is a dimensionless constant ( Frost and Ashby, 1982 ). The effective dislocation core radius r core is assumed to be 5 Å ( Chang and Graham, 1966 ) or 2 b ( Niu et al., 2019 ), and therefore a core = π r 2 core . f l = 1 − f core ≈ 1 , especially when the dislocation density is low. Diffusion coefficients can be expressed empirically as D = D 0 exp (−U/kT ) , where D 0 is a temperature-independent constant, U is an activation energy, k denotes Boltzmann's constant and T represents the absolute temperature in Kelvin. The activation energy for core diffusion U core is approximately 0.6 U l ( Frost and Ashby, 1982 ) and strongly depends on the dislocation character ( Balluffi, 1970;Pun and Mishin, 2009 )  As demonstrated above, at high temperature and low stress, D eff ∼ = D l , and lattice diffusion dominates. While, at higher stresses, or lower temperatures, core diffusion becomes more important. Core diffusion is the dominant mechanism during low temperature creep ( Frost and Ashby, 1982 ) or post-irradiation annealing ( Ferroni et al., 2015 ). Frost and Ashby (1982) demonstrate the range of dominance of this low temperature creep mechanism on the deformation mechanism maps for a wide range of engineering materials. It is the dominant mechanism for a large number of practical applicationsparticularly in the power generating industries. The general formulation in Eq. (1.1) can also be applied to interfaces, free surfaces or grain boundaries, which also provide a short-circuit for the diffusion of point defects and plays a crucial role in diffusional flow ( Foreman and Eshelby, 1962 ) or displacive plasticity ( Li and Ma, 2018 ). Climb-enabled dislocation plasticity has recieved increasing interest in the past two decades. The influence of diffusion on dislocation motion has been considered to investigate loop coarsening and annealing ( Bakó et al., 2011;Clouet, 2006;Liu et al., 2017b;Mordehai et al., 20 08;20 09;Roy and Mordehai, 2017 ), high-temperature creep Danas and Deshpande, 2013;Gao et al., 2017;Haghighat et al., 2013;Keralavarma et al., 2012;Krug and Dunand, 2011;Yang et al., 2015;Zhao et al., 2015 ), interactions at interfaces or free surfaces ( Ayas et al., 2012;Davoudi et al., 2012;Wang et al., 2009 ), and other atypical high-temperature deformation processes Geers et al., 2014;Huang et al., 2014;Keralavarma and Benzerga, 2015;Liu et al., 2017a;Po and Ghoniem, 2014 ). Different numerical methods employed at different time or length scales include, atomistic simulations ( Clouet, 2006;Wang et al., 2009 ), Kinetic Monte Carlo simulations (KMC) ( Kabir et al., 2010;Swinburne et al., 2016 ), discrete dislocation dynamics (DDD) Bakó et al., 2011;Liu et al., 2017b;Mordehai et al., 2008;Thompson et al., 2018;Yang et al., 2015 ), crystal plasticity ( Geers et al., 2014;Petkov et al., 2019 ), Green's function formulations ( Gu et al., 2015 ), and phase-field approaches ( Geslin et al., 2014;. However, despite extensive studies of dislocation climb, only lattice diffusion is accounted for in most of the research, and core diffusion is usually neglected. An understanding of core diffusion processes, including how it influences the motion of dislocations and how this high-diffusivity path contributes to the mechanical behaviour of materials, still remain scarce ( Gao et al., 2011;Niu et al., 2017Niu et al., , 2019Swinburne et al., 2016;Turnbull, 1970;Turunen and Lindroos, 1974 ), even though it is energetically favourable and plays a significant role in the kinetics of diffusion-limited processes, especially at moderate or lower temperature. A systematic investigation of core diffusion and the corresponding self-climb process is therefore essential. As discussed above, a self-climb model is expected to: 1. Contain the essential features of a realistic diffusion process along dislocation cores. 2. Be sufficiently easy to use for the computation of self-climb for complicated dislocation networks. 3. Provide a simple framework to couple with dislocation glide motion allowing simulation of realistic mechanical processes.
Earlier attempts at incorporating core diffusion with dislocation self-climb began with Turnbull (1970) , who investigated the mutual interaction between two prismatic loops. The relative self-climb velocity of the loops was derived analytically and found to be proportional to the driving forces given by Kroupa (1960) . However, dislocation loops were assumed to remain circular throughout the evolution, which is not physical when loops are close to each other or in contact. Turunen and Lindroos (1974) developed a theoretical self-climb model for dislocation lines rather than loops, and proposed that the climb velocity is proportional to the second derivative of the climb force with respect to the dislocation length. Application of this theoretical model is limited due to the difficulty in obtaining the analytical solution of climb forces even for the simplest cases and the model excludes dislocation glide, which usually occurs simultaneously with climb.
As a core diffusion controlled process, self-climb occurs over a time-scale too large to be captured by detailed atomistic simulations ( Zepeda-Ruiz et al., 2017 ) and a characteristic length-scale too small to be described by continuum models. Only recently has self-climb been considered within discrete dislocation (DD) formulations by Gao et al. (2011) . However, this DD model has certain limitations due to some idealised assumptions. For example, vacancy sources/sinks are artificially set for each loop, which is unphysical. The climb velocity is derived to be proportional to the vacancy concentration gradient rather than the chemical potential gradient; consequently the elastic energy release and enthalpy change due to the formation of vacancies are excluded. An equilibrium solution of the vacancy concentration is used in a prescribed uniform vacancy field based on a steady-state diffusion assumption, so that the localisation caused by dislocation curvature is neglected. Furthermore, dislocation glide is not taken into consideration. Only recently, Niu et al. (2017) developed a mesoscopic dislocation dynamics model for vacancy-assisted dislocation climb by upscalings from a stochastic model on the atomistic scale which quantitatively describes the self-climb properties of prismatic loops ( Niu et al., 2019 ).
Solving the full transient diffusion problem to determine the instantaneous climb rate and coupling with dislocation glide is complex and time-consuming, and requires tracking the moving dislocation cores during the evolution of a complex dislocation network ( Liu et al., 2017a ). The variational principle as proposed by Needleman and Rice (1983) can, however, provide a framework to involve multiple thermodynamic and kinetic processes, which allows more efficient numerical schemes to be developed that are suitable for large scale simulations. In fact, a general strategy has been developed by Cocks et al. (1998) and Cocks (1996) to describe the evolution of microstructure due to mass diffusion, including mechanisms that require grain boundary diffusion ( Cocks, 1989;Pan and Cocks, 1995 ) and interface/surface diffusion ( Cocks, 1992;Pan et al., 1997;Suo, 1997 ). The general thermodynamic variational principle is employed in the present work to derive the governing equations for dislocation core diffusion. These are then implemented within the nodal based discrete dislocation dynamics framework to simulate self-climb processes and circumvent the complexity of dealing with the dislocation topology.
The purpose of the present work is then three-fold: 1. Develop a self-climb model to describe conservative dislocation self-climb even at lower temperature. 2. Propose an effective discrete dislocation framework that includes both dislocation glide and climb. 3. Provide a systematic interpretation about the loop coarsening process during post-irradiation annealing.
The paper is organised as follows. Section 2 focuses on the development of a new methodology, to combine the process of core diffusion and dislocation self-climb in a unified framework. The finite element formulation for core diffusion is derived from a variational principle, and then implemented in a discrete dislocation dynamic framework to develop a self-climb DD model. This extends the traditional discrete dislocation framework to allow self-climb motion. In Section 3 , a combined glide/climb DD framework is further developed to enable the simulation of more physically realistic processes. A systematic investigation of loop coarsening behaviour is conducted; including evolution of an isolated prismatic loop, interplay between two loops, and coarsening of a population of loops. The results obtained with this new method are compared with previous experimental and analytic results. Finally, large-scale dislocation patterning is briefly studied.

Methodology
In this section, we develop the governing equations for dislocation core diffusion based on the variational principles proposed by Needleman and Rice (1983) and generalised by Cocks et al. (1998) and Suo (1997) . We then implement these equations into the nodal three-dimensional DD framework to simulate the core diffusion controlled self-climb motion.

The variational principle for core diffusion
The variational principle has been widely used to solve mass diffusion or migration processes in which the rate of evolution of a system is given by the stationary value of the functional. That is, where is a rate potential, which contains contributions from all of the different kinetic and dissipative processes in the system and ˙ G is the rate of change of Gibbs free energy, which is the origin of the thermodynamic driving forces. In this section, we consider a general core diffusion problem along an arbitrarily shaped dislocation line, as schematically illustrated in Fig. 1 , where l is the curvilinear coordinate along the dislocation line, f ( l ) is the Peach-Kohler force acting on the dislocation at position l . The governing quantities symbol and units are given in Table 1 ,  According to Fick's first law ( Shewmon, 2016 ), the constitutive behaviour of all diffusion processes including surface, grain-boundary, and lattice diffusion, follow the same general law; the flux is driven by the chemical potential gradient, where D is an effective diffusivity and μ is the excess chemical potential at position l .
In the problem analysed here, it proves convenient to express the flux as a volumetric rate along the dislocation cores.
For core diffusion along l, μ in Eq. (2.2) then refers to the change in Gibbs free energy per unit volume for the diffusing material along the dislocation. The chemical potential is directly proportional to the climb component of the Peach-Kohler force. This is different from a bulk diffusion-controlled process, where the chemical potential involves contributions from the distribution of vacancies ( Gao and Cocks, 2009 ). The vacancy concentration is not needed in the situations considered in the present paper and the gradient of μ is the driving force for diffusion per unit volume of material, that is, q (l) = −∂μ/ ∂l . Therefore, Eq. (2.2) can be rewritten as, and D ( l ) in Eq. (2.2) refers to the effective diffusivity, where is the atomic volume, D core ( l ) is the core diffusivity, which, in general, is a function of state of the materials, and can be expressed empirically ( Shewmon, 2016 ) as, D core = D 0 exp − U m kT . The pre-exponential factor D 0 is the maximum value of D core (as T → ∞ ), which can be derived based on the physical process of the random-walk problem ( Shewmon, 2016 ). U m is the migration energy for diffusion in the dislocation core, which is found to be strongly depend on the dislocation type (edge, screw, or mixed) by Pun and Mishin (2009) . According to their MD simulations, the rate atoms diffuse along a screw dislocation core is more than an order of magnitude larger than diffusion along an edge dislocation. To reflect this, an effective core diffusion coefficient D core (l) = (D e core sin α) 2 + (D s core cos α) 2 is used in the present work based on the weighted interpolation, with D e core and D s core representing the core diffusion along edge and screw dislocations, respectively.
α is the dislocation character angle between the Burgers vector and the dislocation line direction vector. D core is, therefore, a constant, at a given temperature and dislocation character angle. Meanwhile, the core diffusion process dissipates the energy of the system. The dissipation potential per unit length of dislocation due to core diffusion can be expressed in terms of the diffusion flux j ( l ) along l ( Cocks et al., 1998 ), which is, where the subscript c represents the climb motion. The rate of energy dissipation per unit volume due to dislocation glide, assuming a phonon-drag glide mobility, can be expressed as, where M ( l ) is the mobility for dislocation glide, which is a function of temperature, pressure, and the dislocation type ( Zbib and de la Rubia, 2002 ), v g ( l ) is the glide velocity at position l .
Then the rate potential of the whole system can be derived by integrating ψ( l ) along the total dislocation line length (2.7) By using the Gibbs energy density per unit length, g , the rate of change of Gibbs free energy of the system ˙ G can be expressed as, in which f c ( l ) and f g ( l ) are the components of the Peach-Kohler force f ( l ) at position l in climb and glide directions, respectively ( Gao and Cocks, 2009 ). The Peach-Kohler force is given by, where σ is the total stress tensor, including the externally applied stress and long-range elastic dislocation interactions. The components contributing to climb and glide then read as, with n = (l × b) / | l × b| denoting the unit vector normal to slip plane in the climb direction, m = (n × l) / | n × l| representing the unit vector in the glide direction. Therefore, the ˙ G term can be rewritten as, (2.13) and Eq. (2.1) can be rewritten based on Eqs. (2.7) and (2.13) as, (2.14) It can further be shown that, for any class of assumed compatible flux and velocity fields, there is only one stationary value that minimises the value of . So we can obtain the governing equations for the diffusion process from δ = 0 , which gives the glide and climb velocities of the dislocations.

Mechanism and law for dislocation self-climb dominated by core diffusion
As point defects diffuse along the dislocation core, atoms are added to or removed from a particular part of the dislocation core, resulting in a climb velocity v c which is perpendicular to the slip plane. Conceptually, as illustrated in Fig. 2 , part of the dislocation line climbs by adding a vacancy generating a pair of jogs. Fig. 2 (a) schematically shows the absorption of a vacancy at a pure edge dislocation. Mass conservation requires the volume consumed by climbing, v c bdl , should be equal to the change of the volumetric flux dj , that is, where b is the magnitude of the Burgers vector b . Eq. (2.15) governs the local vacancy flux in the dislocation core when it climbs. The same law works for interstitials but in the opposite climb direction. As stated by Balluffi (1969) , the climb of a mixed dislocation with a substantial edge component should proceed by essentially the same diffusion mechanism as for a pure edge dislocation. Fig. 2 (b) shows the process of the absorption of a vacancy by a mixed dislocation, of which the character angle between the dislocation line direction vector l and Burgers vector b is represented by α. Vacancy absorption occurs with the aid of glissile atomic rearrangements ( Balluffi, 1969;Roitburd, 1965 ), which means that only the edge component of the dislocation line, dl sin α, consumes the absorbed vacancies and the screw component follows by glissile rearrangement. Thus the effective volume consumed during climb should be v c b sin αdl for a mixed dislocation, so that Eq. (2.15) becomes,

dj(l) dl
Note, however, that according to Eq. (2.16) , the climb rate v c → ∞ as the dislocation becomes pure screw, α → 0, which is non-physical. In reality, the glissile rearrangement is expected to be aided by lattice vibrations which are highly temperature-dependent and becomes more difficult as the distance over which rearrangement happens increases. Thus the rearrangements are quite localised at a given range of α. For screw dislocations and dislocations which are close to being pure screw, the glissile rearrangement may not be able to occur. Therefore, to avoid the singularity, a cut-off of the effective dislocation character angle α 0 = 15 • is introduced in the present work. That is, for a dislocation line whose character angle is less than α 0 , the climb rate v c = 0 and motion in the direction of the slip plane normal could be realised by the "pencilglide" mechanism ( Cai and Bulatov, 2004 ). Otherwise, the climb rate v c follows the climb law demonstrated in Eq. (2.16) .

Incorporating the self-climb theory into a 3D-DDD framework
In this section, we implement the variational principle by discretizing the core diffusion formulation into the nodal based 3D-DDD framework to simulate the self-climb process. In discrete dislocation dynamics (DDD), curved dislocation lines are discretized into straight dislocation segments connected by a set of nodes . Each segment has a start and end node. Also, to avoid redundancy, no two nodes can be connected by more than one segment. The approach adopted here for modelling the kinetics of the climb process parallel that described by Pan et al. (1997) for modelling physical processes in which material redistribute by surface diffusion.
Consider a straight dislocation segment (2) bounded by nodes 1 and 2, the degrees of freedom and constraints associated with this segment are illustrated in Fig. 3 . Here, we follow the convention that the nodal values are denoted in uppercase and numbered by superscript, for example, V I c represent the nodal velocity of node I in the climb direction, whereas segment values are indicated in lowercase and numbered by subscript with parentheses, for example b (2) is the Burgers vector of Fig. 3. The degrees of freedom associated with a dislocation segment i . segment (2). J 3 is the diffusive flux across the mid-point of segment (2), denoted as node 3. A local coordinate s is defined. The origin of s is located at the mid-point of the segment (node 3) and the positive direction points from node 1 → 2. s is normalised by half of the segment length l (2) /2, so that s = −1 at node 1, s = 1 at node 2. The unit line direction, ˆ s (2) , and the slip plane normal, ˆ n (2) , for segment (2) are then defined as, (2.17) The nodal climb velocities at nodes 1 and 2, V 1 c and V 2 c , are in the climb directions ˆ n 1 and ˆ n 2 , where we assume that the nodal climb directions are the weighted average of the segment slip plane normals connected to the node, where l (1) and l (2) are the lengths of segments (1) and (2), respectively. Note that, this assumption only applies when the node has only two connections and these two segments have similar orientations. For smooth dislocation curves this is true except for rare cases when the nodal climb direction cannot be explicitly determined. For example, when two adjacent segments are at significantly different orientations or on different slip planes, or when more than two segments meet at a node, in such cases the full (3 component) velocity is used. Assuming a linear variation of the climb velocity along each segment, the climb velocity distribution along segment (2) can be approximately derived as, (2) = ˆ n 2 ·ˆ n (2) . N 1 ( s ) and N 2 ( s ) are the 1D linear shape functions, ( 1 + s ) .
(2.22) According to Eq. (2.16) , along segment (2), mass conservation requires that We define a positive direction for v c when the positive extra half-plane moves upwards or the negative extra half-pane moves downwards. Therefore, a positive climb motion absorbs vacancies or emits interstitials. Also, as mentioned before, only the edge component of a dislocation segment climbs with any connected screw segments following by a quick glissile rearrangement of the atoms ( Balluffi, 1969 ). The climb velocity will be set to zero if the character angle of any segment is less than α 0 . Now we derive the flux j (2) ( s ) by integrating Eq. (2.24) , (2.28) Substituting Eq. (2.27) into Eq. (2.5) and noting that D ( l ) is constant for a given segment, we can derive the contribution of segment (2) to the rate potential ψ c , ] is a 3 × 3 viscosity matrix for core diffusion along segment (2), As illustrated in Eq. (2.35) , the elementary viscosity matrix [ K (2) ] depends not only on the quantities of segment (2), such as l (2) , α (2) and D (2) , but also on how the segment is connected to other segments, through β 1 (2) and β 2 (2) , leading to interaction between elements when assembling the global viscosity matrix.
The total global rate potential c due to climb can be expressed as, is an (i + I) column vector of unknowns including fluxes at the mid-points of all segments and velocities at all nodes, with i denoting the total number of segments and I the total number of nodes.
[ K ] is a global (i + I) × (i + I) viscosity matrix for core diffusion, which is obtained from assembling the elementary viscosity matrices for every segment, The rate of change of Gibbs free energy during climb, ˙ G c , in Eq. (2.13) can also be discretized as, The first term refers to elastic interaction between the (i) segments connected to node I with every other segment ( j ) = ( i ) whereas the second term denotes the elastic self force on node I, due to segments (i); we use the expressions derived by Arsenlis et al. (2007) . The third term is the core force on node I and the final term accounts for an applied stress. The climb force of node I is then, Note that, Eq. (2.23) only guarantees mass conservation in the segment, not at its two ends and hence mass is not conserved at junctions where two or more segments meet. Continuity at each junction node K , requires, where the sum is over all segments ( i ) connected to node K , ˆ l K (i ) is the unit vector along segment ( i ) pointing into or away from the junction node K : ˆ l K (i ) = 1 when j ( i ) flows into the node; ˆ l K (i ) = −1 when j ( i ) flows away from the node. Eq. (2.40) requires that what flows into a junction node is equal to what flow out. A series of Lagrange multipliers, one λ K per junction node K , are therefore introduced to account for the constraints at all junction nodes, so that the functional c can be rewritten as, where the inner summation in the third term is performed over all segments connected to node K and the outer summation is over all junction nodes, K , in the system. Therefore, is a 2 × 3 complementary matrix, where taking segment (2) as an example, (2.45) [ C ] is a I × (i + I) constraint matrix relating the flux j , at the junction nodes, to the velocity and flux degree of freedom [V c ] through mass conservation.
[ λ] = [ λ 1 , λ 2 , . . . , λ I ] T is a I × 1 vector of Langrange multipliers, which physically represent the chemical potentials at the corresponding junction nodes (see Cocks et al. (1998) (2.46) The correct flux and climb velocity pattern is then that which minimises the functional c . So that, δ c = 0 for any small The variational principle in Eq. (2.47) provides the kinetic equations for dislocation climb and the continuity of flux at each node, that is, (2.48) This matrix equation represents a set of linear simultaneous equations which can be solved by using standard matrix methods. The advantage of using Lagrange multipliers to enforce mass continuity at the junctions is that it allows the viscosity matrix [ K ] to be organised into a narrowly banded matrix which allows efficient solvers to be used. Also, the constraints are enforced exactly, that is, mass conservation at the junction nodes are strictly satisfied. We also need to be careful in choosing a solver for Eq. (2.48) , because the leading matrix is not positive definite due to the additional variables introduced by the multipliers. An alternative approach to deal with the constraints without introducing additional variables is through the introduction of penalty function. In the penalty function method, a large positive parameter P is employed in the functional term to penalise solutions which do not conserve mass, and Eq. (2.41) can be written as, (2.49) Theoretically, the constraints are also exactly satisfied as P tends to infinity, which, however, leads to an ill-conditioned system of equations ( Weyler et al., 2012 ) and causes numerical problems. Choosing an appropriate penalty parameter P is, therefore, a balancing act between accuracy, efficiency and stability.
In this paper we concentrate on enforcing the constraints using Lagrange multipliers and determine the rate of climb of the dislocation segments using Eq. (2.48) . Once [V c ] is calculated from Eq. (2.48) , the nodal climb velocity can be derived as, (2.50) The dislocation network can then be updated by V I c t where t is an appropriate time increment. In this process, the lengths of the segments and the connectivity of the nodes will change, leading to complicated topological changes which have been effectively handled in DDD simulations .

Validation and application of the self-climb model
Prismatic loops form after irradiation, enabling a fast core diffusion process. Dislocation self-climb, which is dominated by the core diffusion of point defects, has been proved to be the primary mechanism in the post-irradiation annealing process. For example, Swinburne et al. (2016) demonstrate that the rate of coarsening of an array of dislocation loops of mean radius 25 to 200 Å is around three to four orders of magnitude faster by self-climb at 750K than by vacancy diffusion. Typical diameters of these prismatic loops are up to a few hundred Å and the maximum density is about 10 23 ~10 24 loops / m 3 , which corresponds to a dislocation density of about 10 14 ~10 15 / m 2 ( Kroupa, 1966 ). The coarsening process of prismatic loops after irradiation by a glide and self-climb mechanism is investigated in this paper based on the self-climb model.
Analysis of the movement of prismatic loops prior to coalescence can be divided into three stages: 1. The evolution of isolated prismatic loops. 2. The coalescence of pairs of loops; including parallel and non-parallel loops. 3. Loop coarsening of a large population of loops.

Evolution of an isolated elliptical prismatic loop
For an arbitrarily shaped isolated loop in a crystal, the shape will change to minimise the dislocation energy and balance any forces exerted on it. In this ection, the self-climb model developed above is applied to investigate the evolution of an elliptical prismatic loop. The evolution of the loop would be dominated by core diffusion, in which point defects can only be transmitted along the dislocation core. Therefore, the climb motion only depends on the local atomic rearrangement along the dislocation line and the area enclosed by the loop is conserved.
An interstitial elliptical prismatic loop with an a /2 111 Burgers vector is introduced as an initial configuration, as shown by the blue curve in Fig. 4 (a). The loop has major and minor axes of a 0 = 150 nm and b 0 = 100 nm respectively and is discretized into 36 segments. Parameters for bcc Fe used in the simulation are shown in Table 2 . The core diffusivity along   ( Frost and Ashby, 1982 ) Pre-exponential for core diffusion a core D 0 core = 1 × 10 −23 m 4 / s ( Frost and Ashby, 1982 ) a edge dislocation is, therefore, D e core = D 0 core exp (− U m kT ) . The core diffusivity along a pure screw dislocation is set to be 10 × higher than for an edge dislocation, D s core = 10 D e core , as indicated by available MD simulations ( Pun and Mishin, 2009 ). The Temperature is simulated to be T = 750 K, at which the diffusion process is dominated by core diffusion rather than bulk diffusion  ). An infinite boundary condition is used so that only the elastic dislocation interaction wid etild e F I , and the line tension line (i ) (composed of the elastic self force F I sel f and core force F I core ) are considered as the driving forces (see Eq. (2.38) ).
Snapshots of the loop profile at different times are demonstrated in Fig. 4 (a); a movie of the evolution process is available in the supplementary material S1. The ratio between semi-axes a and b , η = a b , is introduced to describe the shape of the loop. The configuration evolves from an ellipse, with a ratio η = a b = 1 . 5 , to a circle, η = 1 , to reduce the high curvature at the poles and the total dislocation line length. Meanwhile, the configuration converges to a circle with a radius R = 61 . 24 nm = a 0 b 0 , demonstrating the conservation of the loop area and hence conservation of mass. This is also proved by the evolution of the normalised enclosed-area, A / A 0 (with A 0 denoting the enclosed area of the initial ellipse), as the red curve demonstrated in Fig. 4 (b).
An analytic solution for the evolution of η can be derived, by simplifying the ellipse to a rectangular loop with half-height a 0 = 150 nm and half-width b 0 = 100 nm (see Appendix A ), as . Results of the analytic solution for a rectangular loop in Eq. (3.1) and the numerical solution of an elliptical loop given by the DD simulation are demonstrated and compared in Fig. 4 (b). The agreement further validates our self-climb model.

Typical interactions between two prismatic loops
In a multi-loop system, loops can change not only their shapes (the only mechanism possible in a single loop system) but also their relative positions, and even their sizes, to minimise the total energy of the whole system. Typical interactions between two prismatic loops, including parallel and non-parallel loops, are discussed in the present section.

Self-climb of rigid prismatic loops
An earlier attempt to calculate the self-climb coarsening process of two prismatic loops is presented by Turnbull (1970) , which was obtained using an analysis similar to the coalescence and migration of gas bubbles under a temperature gradient ( Gruber, 1967 ). In Turnbull's model, the loops were assumed to remain circular throughout their evolution. A more straightforward derivation of the self-climb motion of a rigid circular loop can also be obtained based on the variational principle developed in Section 2.1 .
Consider the motion of a rigid circular prismatic loop of radius R under the action of forces F x in the climb direction and F z in glide direction, as illustrated in Fig. 5 . The habit plane of the loop is the x − y plane and the Burgers vector is along the z direction. The evolution law for the centers of the rigid prismatic loop in the glide and climb direction can be derived (see Appendix B for details) as, Here, to facilitate our discussion later, we introduce a climb mobility M c for core diffusion controlled climb and a glide mobility M g for viscous glide, so that, where F x and F z can be obtained from their relative position, which will be discussed in the following section. M c and M g are defined as, Noting that the mobility of the loop due to self climb, originally proposed by Turnbull (1970) , was, with D ≡ D core and b 3 , the self-climb mobility M c in Eq. (3.4) agrees well with M c derived from Turnbull (1970) . A similar self-climb mobility was used by Swinburne et al. (2016) to simulate the coarsening of a pair of prismatic loops, where E sc is the activation energy for self climb, a is the lattice constant and ν is the jump frequency. According to the fundmental theory of diffusion ( Shewmon, 2016 ), diffusivity along the dislocation core can be expressed as, Substituting Eq. (3.7) into Eq. (3.6) and using a core 4 a 3 b 2 , gives, Comparing Eq. (3.4) with Eq. (3.8) gives the relationship between the diffusivity D sc used by Swinburne et al. (2016) and the core diffusivity D core used in the present work, To compare with the results in Swinburne et al. (2016) , an effective diffusivity for core diffusion D = D core a core kT = 1 . 0156 × 10 −43 m 7 s −1 J −1 derived from D sc in Swinburne et al. (2016) was used in the following simulations. Now consider the coarsening of two prismatic loops on the same habit plane as shown in Fig. 6 (a). Both analytical calculations, based on the rigid loop assumption, and numerical simulation based on DD were conducted. The initial configuration are the two blue loops with R 1 = 12 nm , R 2 = 20 nm, d = | O 1 O 2 | = 42 nm . Since the elastic force in the glide direction F z = 0 when loops lie in the same habit plane, only self-climb is considered. It is also worth noting that, since the climb forces F x exerted on one loop by the other have the same magnitude, i.e. F 12 x = −F 21 x , the mobility M c is proportional to 1/ R 3 , the self-climb velocity of the loop centres should satisfy v 1 c / v 2 c = R 3 2 /R 3 1 , implying that a smaller loop climbs much faster than a larger one, which is consistent with experimental observations ( Ferroni et al., 2015 ). Other parameters needed for the calculation are shown in Table 2 .
Since the analytic solution Eq. (3.4) was used to calculate the movement of the loop centres, the loops remain circular throughout the simulation. After a simulated time of t rigid = 5 . 358 s the loops touch each other. Snapshots in the x-y plane during loop evolution are shown in Fig. 6 (a), in which the blue line is the initial configuration and the red one is the final configuration. x 1 and x 2 are the displacements of loop 1 and loop 2, respectively. The ratio between the loop While, in the DD simulation, each loop was discretized into 20 segments so that the loops can change their shape. Simulations show that it takes t DD = 2 . 1417 s for the loops to initially touch, which is a faster process than that of the rigid loops. Snapshots of the initial and final configuration of the loops are shown in Fig. 6 (b). The ratio between loop displacements turns out to be η DD = x 1 x 2 = 4 . 2162 0 . 6486 = 6 . 505 , which is larger than η rigid for the rigid loops.
Comparing the final configurations of the rigid and discrete loops in Fig. 6 (a) and (b), we can see that, the discrete loops change their shapes, not only their relative positions, especially when the loops are close, as this lowers the total elastic energy of the system more rapidly, leading to a faster climb rate than for the rigid loops. Meanwhile, local atomic rearrangement is much more energetically favourable than moving the whole loop, especially for a large loop. Therefore, the ratio of the loop displacements for the discretised loops η discrete is larger than for rigid loops η rigid .
The evolution of the total energy is shown in Fig. 6 (d). The black curve is the total dislocation energy E tot of the rigid loops and the blue curve is for the discretised loops. The total energy E tot includes the elastic energy and the core-energy   for an infinite isotropic elastic medium. As shown in Fig. 6 (d), the total energy of the rigid loops decreases slowly until the loops intersect each other, while E tot for the discrete loops keep decreasing and a rapid reduction occurs when the loops coalesce with each other to form a single larger loop. It is found that the equilibrium configuration of these two loops converge to a large circular loop with a radius of R = 23 . 36 nm , as shown in Fig. 6 (c), which satisfies the conservation of loop area with the minimum line length, and the total energy E tot converges to a stable value of E tot = 6 . 71 × 10 −16 J . Kroupa (1960) and Cai et al. (2006) found that the total (elastic + core) energy of a circular prismatic loop is, where r core is the cut-off radius of the dislocation core. Therefore the energy of a loop of R = 23 . 36 nm is E(R ) = 6 . 66 × 10 −16 J , which is also plotted as the red line in Fig. 6 (d) and an excellent agreement is obtained.

Coarsening of two parallel loops by a glide-climb mechanism
Now we simulate the coarsening of two parallel loops in bcc Fe by a combined glide and climb mechanism. The driving force for aggregation is provided by the elastic interactions. Consider two parallel loops of radii R 1 and R 2 with identical Burgers vector, b , as illustrated in Fig. 7 . An analytic expression for the long-range elastic forces exerted, on loop 2 by loop 1, in the glide and climb directions are given by Foreman and Eshelby (1962) , based on the infinitesimal loop approximation, as, where pre-factor Q = Gb 2 / 4(1 − ν) is a constant, d is the distance between the loop centers, α is the angle between the O 1 loop normal and − −− → O 1 O 2 , ranging from 0 to 2 π . w g and w c are geometric coefficients determined by the relative positions of  the two loops, w g (α) = 3 cos α 3 − 30 cos 2 α + 35 cos 4 α w c (α) = 3 sin α 1 + 10 cos 2 α − 35 cos 4 α (3.12) According to Eq. (3.11) , the interaction forces are determined by the relative position of the loops, it is possible to define the extent of the capture zone ( Eyre and Maher, 1971 ), within a given annealing time and local plastic flow stress, to describe the interaction behaviour of loops. Consider a pair of small prismatic loop, with loop 1 at the origin, the force it exerts on loop 2 in the glide and climb directions are schematically demonstrated in Fig. 8 (a) and (b), respectively. For example, when the loops are coaxial, i.e. α = 0 , with loop 2 directly above loop 1, w c (α) = F c = 0 , w g (α) = 24 and F g = 24 Qπ R 2 1 R 2 2 /d 4 > 0 , and the force on loop 1 is equal and opposite ( F g < 0) so the loops would remain coaxial and repel each other by prismatic glide. Whereas when the loops are co-planar, i.e. α = π / 2 , with loop 2 to the right of loop 1, then w g (α) = F g = 0 , w c (α) = 3 , and F c = −3 Qπ R 2 1 R 2 2 /d 4 < 0 (loop 2 would move left and loop 1 would move right). The loops would then remain on the same habit plane and attract each other in the climb direction by self-climb.
However, for the general case when 0 < α < π /2, then both glide and climb mechanisms would operate during the coarsening process. Because the dislocation glide rate is much higher than that of climb, and as the numerical precision is limited, the displacement in the climb direction x c within a time increment would be too small to make any difference to the total displacement, that is x c + x g ∼ = x g . Therefore, bridging the large difference in time scales between glide and climb, and coupling them in a unified computational framework is required to simulate the coarsening process.
As described in the preceding section, the climb velocity can be derived from the finite element formulation in Eq. (2.48) and incorporated into 3D-DDD through the nodal velocity to simulate the core diffusion dominated disloca-tion climb process. The next step is to introduce dislocation glide into this self-climb model. In the present work, a linear phonon-drag mobility law for bcc metals where v J is the nodal velocity of node J, l ( i ) is the length of segment i , with a drag tensor B ( i ) and the summation is over all segments connected to node J . Detailed expressions for the orientation dependent drag coefficient tensor, B , in a bcc material is given in Appendix C . In the present work, pencil glide ( Taylor and Elam, 1926 ) is assumed for the glide motion of the loop, in which glide occurs along 111 directions on any slip plane which contains a slip direction (even non-crystallographic planes). Note that, although a drag coefficient for the climb of a pure edge dislocation B ec is used in B to keep Eq. (3.13) consistent and complete, the contribution to the climb velocity is eliminated by using a very large B ec , with B ec B eg . An effective climb velocity is derived from the core diffusion mechanism from Eq. (2.50) . F J is the nodal mechanical driving force (see Eq. (2.38) ). The elastic interactions between the point defects and the dislocation lines are neglected in the present work.
Once the nodal velocity is computed for every node in both climb and glide directions, the remaining challenge is to bridge the disparate time scale between the two processes. An adaptive time scheme ( Keralavarma and Benzerga, 2015;Keralavarma et al., 2012;Liu et al., 2017a ) is adopted in the present work to perform a simulation over a time duration sufficient for the diffusion process. A small initial time increment t g of ~1ns and only dislocation glide is considered until the plastic strain rate ˙ ε p approaches zero, that is, (i ) is the total area slipped by the n dislocation segments during the time increment t g , and V is the total volume of the computational cell. is a small positive value such that | ˙ ε p | ≤ ensures that a quasi-equilibrium state is reached. Then, the glide process is stopped and climb is activated and a much larger increment t c = 0 . 1 t t g for dislocation climb is used, where t depends on the effective core diffusivity and temperature. It is the time taken for a dislocation to climb one inter-planar spacing distance d (111) = a/ √ 3 . Once the accumulative climb distance of any segment is larger than d (111) , the climb process is stopped and the glide motion is activated and the smaller time increment t g is used again. Both t = 10 t c and t g are set to guarantee that the motion at any time increment is ≈ 5 b . Further details about the adaptive time stepping scheme are given in Liu et al. (2017a) .
We now simulate the loop coarsening process with the climb-glide model outlined above. Dislocation loops in bcc materials are expected to have Burgers vectors of a /2 111 , 100 and perhaps a /2 110 ( Eyre and Bullough, 1965 ). Two interstitial prismatic loops with b = a/ 2 111 and different radii, R 1 = 12 nm, R 2 = 20 nm, were introduced as the initial configuration. The distance between the loop centres O 1 and O 2 was 43.2 nm with a horizontal distance x = 42 nm and a vertical distance z = 10 nm. This is a similar configuration to that analysed in Swinburne et al. (2016) and observed in the TEM by Ferroni et al. (2015) . As shown in Fig. 9 (a), the habit plane of these loops are the x − y plane and the Burgers vector is in the + z direction. The simulated temperature was T = 750 K.
The relative positions of the loops are governed by the elastic interactions in Eq. (3.11) . With an initial angle α 0 = π − arctan ( 42 / 10 ) ≈ 103 . 5 • , the glide and climb forces on dislocation 2 due to dislocation 1 are F g < 0 and F c < 0, respectively, as marked in Fig. 8 (a). The loops would initially repel each other in the glide direction but attract each other in the climb direction. The resulting motion of the loops (loop 2 moves quickly down via glide and slowly left via climb with loop 1 moving in the opposite direction) increases α, until F g = 0 at α 1 ≈ 109.9 • . Glide stops and the loops attract each other by self-climb, further increasing α until α = α 2 when glide occurs again, only now the glide force is attractive, as shown in Fig. 8 (a). Then climb is stopped and the loops glide toward each other decreasing α again until F g = 0 at α 3 = α 1 ≈ 109 . 9 • . This process then repeats, until the loops merge, leaving a larger circular loop. Note that in real materials, Eq. (3.11) would fail to predict the motion of loops when d is so small that loops cannot remain circular (for which no general analytic solution exists), but is a good approximation when the loop separation is much greater than their sizes, d R . The results of a discrete dislocation simulation are shown in Fig. 9 (a)-(d), experimental results for a similar coarsening process reported by Swinburne et al. (2016) is reproduced in Fig. 9 (e).
As shown in Fig. 9 (c), loops can merge with each other even before they move to the same habit plane, because of the non-uniform forces around the loops, neighbouring segments on different loops glide faster than the distant segments, which is not captured with the analytic solution for rigid loops, but is consistent with the experimental observation in Fig. 9 (e). The merged loop in Fig. 9 (c) would then tilt via glide to be on habit plane and allow further climb to form an equilibrium circular loop with R = 23 . 78 nm ∼ = R 2 1 + R 2 2 , as shown in Fig. 9 (d). The displacement of the larger loop is much less than that of the smaller one, especially in the climb direction, which is also consistent with the experimental results. In addition, a recent independent DDD study ( Niu et al., 2019 ) for a similar case showed the time for two circular prismatic loops to coalesce into a single circular loop of R = 94 . 2 b = 23 . 32 nm was 4.1s, which agrees well with the present work. The excellent agreement between our numerical simulation, experimental observations and independent simulations, validates our model and indicates the leading role played by self-climb during the post-irradiation annealing process.

Interaction between two non-parallel loops
Typical interactions between two non-parallel prismatic loops were studied, to validate mass conservation at junction nodes with more than two dislocation segments. The initial configuration is shown in Fig. 10 (a), two non-parallel loops with radii of R 1 = 80 nm and R 2 = 160 nm were introduced with a separation of d = 251 nm on perpendicular habit planes The coalescence process obtained from the DDD simulation is shown in Fig. 10 . The two loops attract each other by self-climb on their habit planes and intersect. Instead of the segments annihilating (as in the case for parallel loops of the same type), the four segments at the intersection form a junction, based on the conservation of Burgers vectors. After which, the evolution of the two loops is driven by the difference in their line tension. The smaller loop shrinks due to its higher curvature ( κ~1/ R ) and hence higher line tension force. Its loss of area is accounted for by growth of the larger loop. A final equilibrium configuration of a circular loop with R = 178 . 76 nm ∼ = R 2 1 + R 2 2 is obtained as shown in Fig. 10 (d), validating that mass conserved is being enforced at the junction node.
In this and the previous sections, we studied simple cases: the self-climb of an isolated loop and pair-wise interactions. For these cases, we can readily obtain analytical solutions for the situation to validate our simulation results. Now we can move on to use our model to simulate patterning within a larger population of dislocation loops.

Coarsening process of a large loop population
The coarsening process of a population of loops is simulated to investigate the emergent interactions and larger-scalepatterning in bcc materials, based on the local rules discussed above. An initial configuration of 60 interstitial-type prismatic loops with identical Burgers vectors, b = a/ 2[111] , were introduced as shown in Fig. 11 (a). These loops were dispersed randomly in an infinite solid over a cubic region of size 20 0 0 b , corresponding to an initial dislocation density of 2 . 8 × 10 14 m −2 . The corresponding loop density is 4 . 99 × 10 20 m −3 , which is similar to experimental data on irradiation-induced loop annealing ( Holmes et al., 1968 ). Radii of the loops range from 15 nm to 50 nm with a random uniform distribution and initial average radius of ~28 nm. Dislocation configurations at different times during the evolution are demonstrated in Fig. 11 , a movie is available in supplementary material S4.
From Fig. 11 , it can be seen that loops repel each other in the vertical direction via glide within a given distance. Whereas loops of the same type aligned along the horizontal direction attract each other via climb and coalesce. Because glide is  much faster than climb, loops repel each other in the glide direction initially, to form a self-organised pattern as shown in Fig. 11 (b). The glide process slows down as the distance between the self-assembled structures increases, and coalescence of loops inside of each self-organised structure is then dominated by self-climb, whereby, loops attract each other via selfclimb to form larger loops. It can also be seen that smaller loops are more mobile than larger ones, which is consistent with Eq. (3.4) .
For closer inspection of the behaviour of the loops, detailed loop interactions, marked as box 1 and box 2 ( Fig. 11 ), are tracked during the evolution. As shown in box 1, closely aligned circular loops in the vertical direction interact with each other by a similar local rule as the coarsening of two parallel loops discussed in 3.2.2 , and form a larger irregular loop. While, horizontally distributed loops, as in red box 2, coarsen into an oblong "finger shaped" loop ( Fig. 11 (b)), and converge to a circular loop over time ( Fig. 11 (c)), which is consistent with experimental observations of the interactions of multiple loops in bcc tungsten by Ferroni et al. (2015) .
Evolution of the total loop area and the average loop diameter are plotted in Fig. 12 . The total area enclosed by the loops, the black curve, is conserved throughout the simulation, with a relative error smaller than 0.13%, which is caused by topological remeshing of the dislocation network. As the coarsening process continues, the number of loops decreases, and as the total loop area is constant, the average loop size increases over time; shown by the red curve in Fig. 12 . This agrees qualitatively with the experimental data from Ferroni et al. (2015) . It is also worth noting that the rate of increase of the average loop size is much larger during the early stage when glide motion dominates. Once the dislocation cell has F. Liu, A.C.F. Cocks and E. Tarleton / Journal of the Mechanics and Physics of Solids 135 (2020) 103783 formed, after ~10s of annealing, the distance separating cells becomes larger and the coarsening process inside of each cell is dominated by the self-climb motion, which is much slower.

Concluding remarks
Since the mean jump frequency of atoms in the dislocation core is much higher than for atoms in the regular lattice, the diffusivity along the dislocation core is much faster. This high diffusivity path, provided by dislocations, plays an essential role in the kinetics of diffusion-limited processes; especially at or below half the melting temperature ( Frost and Ashby, 1982 ). Investigations of this process remains scarce, particularly in relation to post-irradiation annealing.
In the present work, a new methodology has been proposed for describing the core diffusion dominated dislocation self-climb process. Major contributions are 1. Deriving the governing equations for core diffusion from a variational principle for the evolution of microstructures, as demonstrated in Eq. (2.48) . 2. Guaranteeing flux conservation during dislocation reactions such as during junction formation, by introducing a series of Lagrange multipliers. 3. Discretizing the continuum core diffusion formulation for implementation into the nodal-based 3D-DDD framework to simulate the self-climb motion; which allows traditional discrete dislocation modelling to be extended to a new class of physical problem.
The evolution of an isolated loop and the interactions between a pair of parallel loops are used to validate the method by comparison with analytic solutions and experiments.
This glide/climb model, was used to study the interactions and large-scale patterning in bcc materials, during postirradiation annealing by simulating loop coarsening. Results demonstrate that during the early stage, loops repel each other in the glide direction to form self-assembled dislocation cells. Then, loops inside of the cells attract each other via self-climb and coarsen into larger loops over a much longer time scale.
Despite the idealised assumptions, for example, the cut-off angle for the climb of near screw dislocations, or the decoupling of the glide and climb in each increment, the proposed methodology sheds new light on the simulation of core diffusion dominated mechanical processes and enables physically based predictions of post-irradiation annealing with high computational efficiency. This new method will be used in future work to investigate low-temperature creep in single crystals ( Frost and Ashby, 1982 ). Another interesting topic is the so-called displacive plasticity for micro-and nano-sized crystalline materials ( Li and Ma, 2018 ), where the strength is controlled by diffusional plasticity. We would like also to emphasise that although the present study focuses on the post-irradiation annealing process in a bcc Fe system, the underlying principle, and the self-climb/glide mechanism, are general and could be used to simulate dislocation activity and inelastic deformation in many different materials.

Declaration of Competing Interest
None.
where F is the Peach-Koehler force per unit length in the climb direction. In addition, mass conservation requires j a (y ) + j b (x ) = 0 , so that y ˙ ( where A = 3 DF / (a 0 b 0 b 2 e ) . Eq. (A.5) provides the evolution of the ratio η = y/x, with η → 1 for large t , resulting in a square loop.

Appendix B. Mobility law for a rigid prismatic loop
Consider a rigid prismatic loop of radius R under the action of forces F x in the climb direction and F z in glide direction, as illustrated in Fig Due to symmetry, j(θ = 0) = 0 at s = 0 therefore the constant of integration is zero. Now, the variational functional of the whole loop can be derived by substituting Eq. (B.6) into Eq. (2.14) . Assuming that the effective diffusivity D and the viscous coefficient B are constant at a given temperature, we have where l is the unit vector in the line direction, I is the identity tensor. B l is the drag coefficient along the dislocation line direction, with B l B s to allow nodes to move along the line direction freely. The drag coefficient tensor of a non-screw dislocation segment is calculated based on an interpolation function ( Cai and Bulatov, 2004 ), where B eg and B ec are the drag coefficients for the motion of a pure edge dislocation in the glide and climb directions, respectively. m ≡ n × l is the unit vector in the glide direction, with n denoting the unit vector normal to the slip plane.
In the present work, B eg < B s due to the fact that the mobility of a pure edge dislocation is higher than that for a pure screw. The contribution to the climb velocity is eliminated by using a very large B ec , with B ec B eg . The actual climb velocity is derived from the core diffusion mechanism, Eq. (2.50) . The temperature dependence of the drag coefficients ( Po et al., 2016;Queyreau et al., 2011 ) is neglected in the present work. Which is an appropriate simplification considering the huge difference between glide and self-climb rates. B l = 10 −4 B s and B ec = 10 10 B eg are used to satisfy the above constraints. B eg = 5 × 10 −4 Pa · s and B s = 1 × 10 −2 Pa · s are adopted in the present work from the atomic simulations conducted by Wang and Beyerlein (2011) .