An improved method to model dislocation self-climb

Dislocations can provide short circuit diffusion paths for atoms resulting in a dislocation climb motion referred to as self-climb. A variational principle is presented for the analysis of problems in which fast dislocation core diffusion is the dominant mechanism for material redistribution. The linear element based self-climb model, developed in our previous work [1] Liu, Cocks and Tarleton (2020 J. Mech. Phys. Solids 135 103783), is significantly accelerated here, by employing a new finite element discretisation method. The speed-up in computation enables us to use the self-climb model as an effective numerical technique to simulate emergent dislocation behaviour involving both self-climb and glide. The formation of prismatic loops from the break-up of different types of edge dislocation dipoles are investigated based on this new method. We demonstrate that edge dipoles sequentially pinch-off prismatic loops, rather than spontaneously breaking-up into a string of loops, to rapidly decrease the total dislocation energy.

(Some gures may appear in colour only in the online journal)

Introduction
The diffusion of atoms in crystalline materials is a thermally activated process which is drastically enhanced by lattice defects [2], such as surfaces, grain boundaries [3], and dislocations [4,5]. The short-circuit paths provided by these defects are thought to originate from the loose-packed and highly disordered atoms that lower the activation energy for diffusion and enhance the atomic jump frequency. When carried by dislocations, the diffusion mechanism is known as core diffusion or pipe diffusion, as the core of dislocations act as pipes for atoms to diffuse down. However, while surface diffusion and grain boundary diffusion have been studied and discussed in detail [6][7][8][9][10][11][12], only limited results are available for pipe diffusion along dislocation networks [13]. It is challenging to extract core diffusivity from experimental observations [14], since core diffusion often occurs simultaneously with lattice diffusion and occupies far fewer atomic sites. Core diffusion is, therefore, usually neglected in most studies [15][16][17][18] or treated phenomenologically as an enhancement to bulk diffusion [19]. However, at lower temperatures, core diffusion plays a dominant role in mass diffusion, particularly at high stresses. The recent direct observation of core diffusion effects [5,20,21] demonstrated that the core diffusivity in aluminium can be increased by three orders of magnitude near 600 K [5], and the difference can be up to six orders of magnitude in bcc Fe at 700 K [4]. The rapid atomic diffusion along dislocations represents in itself an interesting physical phenomenon and plays a signi cant role in a wide variety of material behaviours, such as creep [19,22], loop or particle coarsening [4,[23][24][25][26][27], precipitation and phase transformations [5,28,29], dynamic strain ageing [30,31], and solute segregation [32][33][34]. However, the role of core diffusion on dislocation motion has received little attention and is not well understood.
During the core diffusion process, atoms can be absorbed or emitted by the core region, leading to climb out of the original slip plane, known as dislocation self-climb. The selfclimb of dislocation loops, as rst proposed by [35], differs from the commonly studied non-conservative climb [18,[36][37][38][39] in the sense that atoms can be transferred only along the perimeter of the loop and the total area of the loop projected perpendicular to its Burgers vector remains unchanged, namely, the relaxation volume of loops stay constant [40]. In self-climb, atoms are rapidly rearranged in the dislocation core region to balance the difference of chemical potential along the dislocation line, which is more energetically favorable than mass exchange with the surrounding lattice (bulk diffusion), particularly at low temperatures where bulk diffusion could be many orders of magnitude slower than core diffusion [4]. Also, the analysis of core diffusion-controlled mechanisms tends to be more straightforward; for example, the three-dimensional lattice diffusion problem [15,17,36,38,39,41] can be reduced to the consideration of one-dimensional diffusion along the dislocation line. Earlier attempts at analysing dislocation self-climb mainly focus on the transport of prismatic loops, due to the simplicity in dislocation character, and their importance in irradiated materials. Following the pioneering theoretical studies on the dislocation self-climb of prismatic loops in the 1960s [35,42,43], Turnbull [26] constructed a model for loop motion by self-climb, in which loops remain circular throughout the process. This theoretical analysis is further supported by the recent experimental observation [4] that the self-climb velocity of a circular prismatic loop is proportional to the driving force and inverse cube of the loop size; indicating higher mobility for smaller loops as observed in post-irradiation annealing experiments [44].
A simple derivation for the transport of a rigid prismatic loop by self-climb was proposed in our recent study, based on a variational principle [1]. The theoretical analyses of self-climb mobility, with the rigid loop assumption, provide a description of the coarsening process of prismatic loops, it fails when loops are close or in contact. In real materials, loops change both their positions and shapes to minimise the total dislocation energy during the coarsening process while conserving the total enclosed area [1,45]. More importantly, the transport of a prismatic loop on its habit plane is just a special case of self-climb of a pure edge dislocation. For the motion of an arbitrary shaped dislocation line, a more general self-climb model is needed to couple atomic diffusion with dislocation self-climb motion.
Molecular dynamic (MD) simulations are well suited to study core diffusion with atomic resolution [46]. MD has been used to calculate various core diffusion dominated behaviours including: the core diffusivity along different dislocation types [14,47], the annihilation of prismatic half loops via pipe diffusion during nanoindentation [48], and the conservative climb motion of a cluster of self-interstitial atoms towards an edge dislocation [49]. Collective dislocation motion involves both climb and glide, making it far beyond the length and time scales accessible with atomistic simulations. A solution is to develop mesoscopic coarse-grained models that operate at intermediate length and time scales, providing a bridge between atomistic and macroscopic models. Geslin et al [50] proposed a multiscale approach to model dislocation climb of jogged dislocations with consideration of both bulk and core diffusion between jogs. The analytical climb rate was derived by assuming a periodic distribution of jogs along the dislocation line. The climb rate was then upscaled to a phase eld model to simulating dislocation climb at a larger scale. Later, Niu and co-authors [45,51] developed a discrete dislocation dynamics (DDD)-based climb model by upscaling from a stochastic jog dynamic model on the atomic scale. Both the pipe diffusion equation with Dirichlet boundary conditions at the jogs and the bulk diffusion equation with Robin boundary conditions near the dislocation were solved in [51], to derive the jog mobility. Note that, among these multiscale models, dislocation climb motion is represented by the dynamics of the pre-existing jogs (neglecting nucleation), which act as perfect sinks or sources of vacancies. This is a good assumption for bulk diffusion when the jog density is high enough to maintain an equilibrium vacancy concentration along the dislocation core [17]. But for a low jog density, as is likely at lower temperatures and higher stresses, the stresses tend to sweep the jogs away and push them towards the dislocation ends. Consequently in this regime the nucleation of jog-pairs needs to be taken into consideration [52].
In our recent work [1], a self-climb model was proposed, which employs a variational principle for the evolution of microstructure [53,54]. In this model, a nite element based analysis for the one-dimensional core diffusion process was implemented. A dislocation self-climb model was then developed by implementing this core diffusion formulation into the nodal based three-dimensional discrete dislocation dynamics (DDD) framework [55], to extend the traditional DDD method to simulate self-climb and glide of any arbitrarily shaped dislocation network.
However, the DDD method is still limited by its high computational cost [56]. Improving the computational ef ciency [57,58] is, therefore, a challenging but necessary task. In a general DDD simulation [59,60], the most time-consuming parts include: (1) the time integration of the equation of motion, which has been effectively sped up ∼ 100× by time subcycling [57,58]; (2) the calculation of seg-seg interactions between different dislocation segments, which is an N-body problem and increases dramatically with strain. This problem has been accelerated by a parallelization strategy on the graphical pgrocessing unit (GPU) [58,61]. In the newly developed DDD-based self-climb model [1], the additional time-consuming part is the calculation of nodal climb velocities. In this method, dislocations are discretised into a series of straight segments (1D nite elements). The climb velocity is de ned at each node and varies linearly along the segment, and the diffusive ux is de ned across the mid-point of each segment. To enforce the ux continuity at nodes where two segments meet, a series of Lagrange multiplier are introduced, one for each node. A set of linear simultaneous equations is then derived as the kinetic equations for self-climb [1], is the (N + m) × 1 zero vector. The advantage of using Lagrange Multipliers to enforce the ux continuity at junctions is that it allows the [K ] matrix to be narrowly banded and effectively solved with standard matrix methods. However, in (1) the leading matrix of the equations is then not positive de nite, due to the introduction of Lagrange multipliers, leading to a signi cant increase in the computational cost which becomes problematic as the number of segments increases. Consequently the purpose of this work is to improve the computational ef ciency of the self-climb model, to enable the simulation of collective dislocation evolution involving both self-climb and glide motion.
The present paper is organised as follows. In section 2, a new discretisation method is developed by the adoption of the paired-linear element based method (PLEBM), to circumvent the usage of Lagrange Multipliers for every node. The PLEBM is validated and compared with the linear element based method (LEBM) to show the improvement in capability and efciency. Finally, in section 3 new physical insight into the break-off of dislocation dipoles, where core diffusion-controlled self-climb is the dominant mechanism, are obtained based on the PLEBM.

Methodology
The variational principle for core diffusion developed in our previous work [1] is brie y reviewed here. We then implement this formulation with a paired-linear element to derive the kinetic equations for self-climb.

The variational principle for core diffusion
Variational principles have been rigorously formulated for a series of mass transport mechanisms, including grain-boundary diffusion [62], surface diffusion [63], grain-boundary migration [64] and coupled grain-boundary and surface diffusion [65]. The essential idea carries over to dislocation core diffusion. That is, among all the virtual velocities of microstructures and virtual diffusive uxes that satisfy matter conservation, the actual velocity and ux elds minimize the functional of the system Π [65], where Ψ is a rate potential term involving contributions from all possible dissipative processes. G is the energy rate term, which is the origin of the generalised thermodynamic force. In (2) it is the combination of thermodynamics and kinetics that determines the actual evolution path and the nal state. Here, the one-dimensional core diffusion-controlled self-climb process is investigated.
The general situation analysed in the current paper is schematically shown in gure 1. An arbitrarily shaped three-dimensional dislocation network, in an in nite domain, is subjected to long-range elastic interactions between dislocation segments. Dislocations can either glide in their original slip plane, or relocate themselves by the diffusional transport of material to move perpendicular to their original slip plane. As stated in the introduction, we consider the situation of lower temperatures when lattice diffusion is negligible in comparison with core diffusion; such that mass can only diffuse along the dislocation line in the core area.
Here, we state two basics laws before developing the variational functional of the system: The law for mass diffusion. We adopt a classical description of the diffusion process, the volumetric ux j(l ) driven by the gradient in the chemical potential as described by Fick's law, where l is curvilinear coordinate along the dislocation line, D core is the core diffusivity, a core is the cross sectional area of the dislocation core, k is the Boltzmann's constant, T is the absolute temperature and µ is the excess chemical potential; For core diffusion along dislocation lines, µ in (3) then refers to the change in Gibbs free energy per unit volume for the diffusing material along the dislocation; it is directly proportional to the climb component of the Peach-Koehler force. The gradient of µ drives the diffusion of material with a unit volume.
is the effective diffusivity. Along the dislocation line, both the local stress and concentration of point defects contribute to the chemical potential µ.
The law for dislocation climb. As mass diffuses along the dislocation core, atoms may be deposited or removed from the extra-atomic-plane, leading to climb motion perpendicular to the slip plane. During this process, mass conservation requires the climb velocity is balanced by the variation of the ux, where dl sin θ is the edge component of the dislocation line increment dl, θ is the angle between the Burgers vector b and dislocation line vector l. As discussed in [1], this ensures only the edge component of a dislocation line can climb and consume the atomic ux. We now develop the variational functional of the system, as derived in [1], where f c is the Peach-Koehler force component in the climb direction, v c is the climb velocity and L is the total dislocation line length in the system. The rst term on the right-hand side of (6) is the rate potential due to core diffusion, and the second term is the rate of change of Gibbs free energy during the self-climb motion. Different virtual motions give different values of Π. Of all virtual motions, the actual motion renders Π stationary, that is,

Implementation of the variational principle
The kinetics, which here is atomic diffusion, proceeds locally. The idea is to transcribe local information onto the total energy landscape, to nd the global equilibrium state. The continuum structure is approximated by discrete elements and a viscosity matrix assigned to every point on the total landscape. In using the variational principle, it is impractical to search for the exact motion from all virtual motions. Rather, the solution is selected from a restricted set of virtual motions, and the solution which minimises the variational functional Π approximates the exact motion.
We have described how numerical schemes can be developed from (6) and (7) for threedimensional dislocation networks in our previous work [1]. We limit our attention here to the differences in the paired-linear element discretisation compared to the previous linear element [1]. As demonstrated in gure 1, dislocation lines are discretised into a series of straight segments (1D elements). Typical types of dislocation connections are shown in gures 1(a) and (b).
We initially restrict our attention to the majority of situations, demonstrated in gure 1(a), where each node connects to only two segments. We will see in the following discussion that it proves convenient to isolate a segment, for example, segment (13) in gure 1(a). Along the segment, the nodal climb velocity V c and the diffusive ux J are de ned at each node. Here, we follow the convention that the nodal values are denoted in uppercase and numbered by subscripts, whereas segment values are indicated in lowercase and numbered by superscripts with parentheses. For example, V c1 is the nodal climb velocity at node 1, and J (12) 1 is the ux at node 1 along segment (12). b (12) , l (12) and n (12) denote the Burgers vector, unit line vector and unit normal vector of segment (12), respectively. In the following discussion, the functional Π is discretised in terms of the degrees of freedom de ned above.
To avoid using Lagrange multipliers at every node, and avoid the non-positive coef cient matrix in the kinematic equation (1), a different discretisation method is used here. An additional dummy node is located at the midpoint of each segment, node 2 in gure 1(a). To standardise the following computation, a local coordinate s is de ned. The origin of s is located at the start point of a given segment, and the positive direction points from node 1 → 2. s is normalised by the segment length l (12) , so that s = 0 at node 1 and s = 1 at node 2. Assuming a linear variation of the climb velocity along both segment (12) and segment (23), so that the climb velocities along segments (12) and (23) are given respectively as, where s ∈ [0, 1], β 1 = N 1 · n (12) , β 2 = N 2 · n (12) = N 2 · n (23) , β 3 = N 3 · n (23) and N 1 , N 2 , N 3 are the directions of the climb velocity at nodes 1,2 and 3, respectively. These are assumed to be the weighted average of the segment slip plane normals connected to the node, for example, (12) n (12) + l (23) n (23) |l (12) n (12) + l (23) n (23) | .
Not to be confused with N (12) 1 (s) and N (12) 2 (s), the isoparametric shape functions for segment (12), N (23) 2 (s) and N (23) 3 (s) refer to segment (23), Mass conservation requires, and so, and the ux is j (23) where l e = lsinα is the edge component of the segment length l. Note that there are now 7 degrees of freedom for segment (13), Flux continuity at node 2 requires that the ux into node 2 must ow out, so that J (12) 2 = J (23) 2 = J 2 . This removes one of the degrees of freedom. Flux continuity is required at every node. Therefore, for a node connected with two segments, only one degree of freedom is needed for the ux. One degree of freedom, J I , is therefore used for the ux at each node I. Care is required for triple or quadruple junction nodes with more than 2 segments, such as node 4 in gure 1(b). In such cases, additional degrees of freedom are needed for the ux along different segments. This is discussed and evaluated further in section 2.4.
[K ] is the global viscosity matrix assembled from all elementary viscosity matrices [k (i) ].
The rate of change of Gibbs free energy during climb motion can also be discretised as, where N is the total number of dislocation nodes, V cI and F cI are the nodal climb velocity and nodal climb force at node I. [F c ] is an N × 1 nodal climb force vector corresponding to the global nodal velocity vector [V c ]. As stated in our previous work [1], the nodal climb force on node I, F cI , is where F I is the full nodal force at node I, with normal N I , as de ned in (10).F I is the elastic interaction between the (i) segments connected to node I with every other segment (j) = (i). F self I refers to the elastic self force on node I, due to segments (i). F core I is the core force on node I and F app I accounts for an applied stress. The variational functional Π can now be expressed in terms of the discretised degrees of freedom by substituting (22) and (23) into (6), Substituting (25) into (6), are the kinematic equations for dislocation self climb motion. A set of linear simultaneous equations with a narrowly banded matrix [K ], which can be solved with standard matrix methods. Once (26) is solved, the generalized coordinates can be updated for a small time increment. The process is repeated for many increments to evolve the dislocation con guration.

Validation of the PLEBM
Importantly, (26) is preferable to (1), derived with the linear element based method [1], in terms of improving the coef cient matrix by eliminating the zero diagonal elements. We now compare the accuracy and computational cost of the paired-linear vs linear element. A benchmark study was performed of the annealing of an isolated elliptical prismatic dislocation loop in an in nite domain in α-iron, at T = 750 K. The initial semi-major axis is a 0 = 400 nm and the semi-minor axis is b 0 = 100 nm (the initial aspect ratio a 0 /b 0 = 4). Other parameters of α-iron are obtained from table 2 in [1]. The loop is discretised into 30 straight segments. The evolution of the loop's aspect ratio a/b over time is shown in gure 2(a), the black curve corresponds to the result obtained with the PLEBM, and the red curve with the LEBM. The corresponding loop pro les during some points in the evolution are also shown. The loop converges to a circle with R ≃ 200 nm to conserve the enclosed loop area. The accuracy of the PLEBM is demonstrated by the excellent agreement between the two methods as shown in gure 2(a), and the consistency with the previous theoretical [42] and numerical [45] studies.
It also worth noting that, although the kinematic equations in the LEBM (1) and the PLEBM (25), give the same results, there is a signi cant reduction in the computational cost achieved with the PLEBM; as shown in gure 2(b). The computational time required to calculate the climb velocity increases with the segment number with both methods. However, the growth rate of the LEBM is much higher than the PLEBM, which becomes signi cant as the segment number increases.

Junction nodes
In the derivation of the kinematic equations for climb velocity, we have implicitly assumed that: (1) every node only connect with two segments and (2) every segment connected to the same node are on the same habit plane, as in gure 1(a). However, for triple or quadruple nodes with more than two connections, as in gure 1(b), (25) and (26) are no longer valid. At a junction node, the direction of the climb velocity N is unknown and (10) is no longer applicable. In this case, a full nodal velocity with 3 components must be used, which is required for nodes connected to segments which are on different habit planes. Moreover, ux continuity at a junction node requires mass owing into the junction node must ow out. We therefore assign more degrees of freedom for the nodal ux, and further constrain the range of admissible elds by requiring that the ux of material into a triple or quadruple junction node is zero; i.e. at each junction node where the summation is over every segments i connected to node I. m (i) = ±1 so that, m (i) = 1 when J (i) ows away from node I, and m (i) = −1 when J (i) ows into node I. We account for these constraints by employing a series of Lagrange multipliers, one per junction node, so that the functional Π in (25) can be extended to include junction nodes as, where the outer summation in the last term is performed over all triple and higher order junction nodes, I, in the dislocation network. Taking node 4 in gure 1(b) as an example, the inner summation can be expanded as, So that, where [λ] is a P × 1 vector of Lagrange multipliers, which are interpreted as the chemical potential at junction nodes, with P indicating the total number of triple and quadruple junctions and [C s ] is a complimentary matrix. We now reconstruct the functional Π as, Taking the stationary value of Π, for any small perturbations Thus (26) can be rewritten to derive the general kinematic equations, Although the structure of (34) is similar to (1) from the LEBM, the size of the [0] matrix in the leading coef cient matrix is now much smaller than in (1). Because multipliers are only introduced for the triple and quadruple junction nodes, which are a small proportion of the total number of nodes. Once (34) is solved, we can then obtain the nodal climb velocity and update the dislocation network for a small time increment, and remesh the topological connections [66] based on the nodal DDD framework.

Applications
In this section, we describe how the variational principle presented in the previous section can be used to simulate dislocation evolution in engineering materials. The examples described in the following subsections are chosen to illustrate the proposed method, and to demonstrate that it produces accurate results by comparison with analytical results and available experimental observations reference [13].

Analytical solution for self climb of elliptical prismatic loop
It is of interest to construct a simple analytical model for the evolution of an elliptical prismatic loop, rather than the rectangular loop derived previously [1]. Sketched in gure 3 is an elliptical prismatic loop in an in nite elastic domain. The two semi-axes satisfy a 0 b 0 . For a dislocation core diffusion problem, the enclosed area of the loop is conserved as the shape of the loop changes. Assume that the loop remains an ellipse during the evolution. Let r 0 be the radius of a circle having the same volume as the elliptical loop. This process can be described by, where α is the shape parameter, with a/r 0 = r 0 /b = α 1, which evolves with time; α = 1 corresponds to a circular loop and α → ∞ corresponds to an edge dislocation dipole. At a xed time, (35) and (36) trace the entire ellipse as θ varies in the interval [0, 2π]. Let r be the position vector of a point (x, y) on the loop, which has an outward unit vector normal n. As atoms diffuse along the prismatic loop in the core region, the dislocation segment moves at a velocity v c = n ·ṙ, whereṙ = dr/dt. With r = (x, y), r = (r 0 α cos θ, r 0 α −1 sin θ) (38) r = (r 0α cos θ, −r 0 α −2α sin θ) So that (37) can be rewritten as, Denote j s as the volumetric ux of atoms due to core diffusion. Mass conservation requires that, where b e is the edge component of the Burgers vector b, for a prismatic loop, b e = b and where the integration constant C = 0 due to the symmetry condition, j s (θ = 0) = 0. Then the rate potential of the whole loop Ψ becomes, where L = 4aE(k) is the circumference of the ellipse, with E(k) denoting the second complete elliptic integral with modulus k = √ a 2 − b 2 /a. The complete elliptic integral of the rst and second kind, K and E, with modulus k = √ a 2 − b 2 /a are, So that, De ne, So that Ψ can be expressed as, while the energy of an elliptical prismatic loop [67,68] is given as, where r c is the core radius of the dislocation. The rate of change of Gibbs free energyĠ is, where I 3 = dE el /dα is given in B. The variational functional now becomes, Taking the stationary value, which gives,α With the initial value α 0 = a 0 /b 0 , the ODE, (56), gives the evolution of the shape parameter α, to describe the loop pro le. Changes in total dislocation energy tend to reduce the length and curvature of the loop, and mass conservation maintains the area enclosed by the loop, driving the elliptical loop into a circular pro le. Next, we compare the results obtained from the analytic solution in (56) with the results from the PLEBM by conducting DDD simulations. Snapshots of an isolated elliptical loop with different initial aspect ratios: α 2 0 = a 0 /b 0 = 4 and α 2 0 = a 0 /b 0 = 16, evolving from an ellipse into a circle, are shown in gure 4. It is expected that the loop will evolve into a circular loop with radius √ a 0 b 0 if mass is conserved, which is found to be the case with both the analytical and DDD results.
With a 0 /b 0 = 4, the analytic solution agrees well with the DDD solution. However, as the initial aspect ratio increases to a 0 /b 0 = 16, with DDD the ellipse initially evolves into a boneshape to reduce the high curvature at the tips of the ellipse, and converges into a circle gradually, which is much more energetically favorable than the process described by the analytic solution. This is because, in the deviation of the analytic solution, we implicitly assume that the loop remain an ellipse, so that a/r 0 = r 0 /b during the evolution. This is a good approximation for the loop evolution to some extent, which was also employed by Sun and co-authors [69] and Cocks [70] to describe the growth of an ellipsoidal void controlled by surface diffusion. However, it is too constrained to capture the evolution of an ellipse with high aspect ratio. The additional degrees of freedom in the DDD method allows it to more accurately simulate the physical process. This raises the question: how do loops evolve if the aspect ratio increases further?

Formation of prismatic loops from dislocation dipoles
Dislocation dipoles and multi-poles are ubiquitous features of the dislocation debris in deformed crystals. Edge dipoles are commonly reported [13,71,72] as arising from the dragging out of jogs on screw dislocations, called jog dragging, as sketched in gure 5(a), or by trapping edge dislocations of opposite sign moving on parallel basal slip planes, called edge trapping, as shown in gure 5(b). Theses dipoles are sessile, they remain as prominent debris in deformed crystals, because they cannot self-annihilate by glide. Consequently they play a signi cant role in strain hardening and patterning in stressed metals. At high temperature or during annealing of deformed crystal, these dipoles can annihilate and break up into prismatic loops by a conservative dislocation climb mechanism, which is primarily controlled by core diffusion along the dislocation line. In certain cases, strings of prismatic loops can form [73]. Attempts to study the annihilation of the dipoles mainly focus on an energetic analysis, in which the dipole uctuation energetics are evaluated as a function of uctuation wavelength and amplitude with either a square wave [74] or sinusoidal wave uctuation [71]. Detailed loop formation mechanisms, by spontaneous break-up or sequential pinch-off from the dipole, still remain controversial [71,74,75]. It is therefore of interest to simulate the formation of prismatic loops from the break-up of dislocation dipoles or elongated loops with DDD.
In the following subsection, we examine two simple representative situations as schematically shown in gure 5, with the mechanism of dipole annihilation of edge dipoles simulated using the PLEBM. We limit our attention initially to the pure self climb process of jog dragging dipoles shown in gure 5(a). We then consider the breakup of edge trapping dipoles shown in gure 5(b), which involves both glide and climb.
3.2.1. Formation of prismatic loops due to pinch-off from the end of a closed-end edge dipole.
A jog dragging dipole as shown in gure 5(a) is introduced as the initial con guration. The typical height of dipoles observed are 10 to 100 nm. The separation distance between two dislocations making up the dipole is therefore set at H = 20nm and the length of the dipole is set at L = 800 nm. The Burgers vector b = 1/2 [111]. The evolution of the total dislocation energy is plotted in gure 6(a) as a function of time, and the snapshots of loop pro les at points marked A-D are demonstrated in gure 6(b); available experimental observations from [13] are presented in gure 6(c); a movie of the evolution process is available in the supplementary material S1.
Our simulation results show that loops are pinched-off sequentially from the tip of the dipole rather than the dipole spontaneously breaking-up into a string of loops. This is consistent with the available experimental observations [13,73]. Initially, the diffusion distance is short, as material ows much faster at the tips than at the middle of the dipole due to the higher curvature. The near tip region thus contracts due to mass conservation as shown in illustration B in gure 6(b). This leads to the rst pinch-off from the two ends. The process repeats itself until the aspect ratio of the remaining loop in is too small to pinch-off, which will converge to a circle as described in gure 4. The evolution of the dipole is a balance between the dislocation self-energy and the elastic interaction energy. Forming loops would decrease the dislocation line curvature which decreases the elastic interaction energy, but increase the total line length which reduces the dislocation self-energy. It is also worth noting that each pinch-off event is accompanied by a signi cant decrease in energy, indicating that the pinch-off of loops is energetically favorable. It is noteworthy that the evolution of this initial con guration, as shown in gure 7(a), involves both glide in the slip plane and climb perpendicular to the original slip plane. A linear phonon-drag mobility law for bcc metals is, therefore, employed to calculate the glide velocity. To bridge the large time scale separation between glide and climb motion, as stated in the introduction, an adaptive time scheme [1] is adopted here. Thus, different time scales are used, small time increments ∼ ns are used to deal with the glide steps, whereas much larger time increments ∼ ms are used for the climb steps. Figure 7 shows snapshots of the dislocation evolution; a movie of the evolution process is available in the supplementary material S2. There is approximately six orders of magnitude difference in the time between the different stages, which signi es the dominant mechanism (glide or climb) at different stages. Initially, the dislocations glide towards each other on their original slip planes at a very fast rate, until they are trapped by each other, as shown in gure 7(b). This mechanism of dipole formation seems to be very frequent at the beginning of plastic deformation when a large number of dislocation loops expand on parallel slip planes. This glide process only takes 9µs and results in a stable edge dipole. Then, the attraction from the anti-parallel edge segments facilitates self-climb normal to the original slip planes, at a much slower rate, as shown in gure 7(c). This climb motion gives rise to the uctuation of the dislocation line to conserve mass. In addition, according to (5), a mixed dislocation climbs a larger distance when consuming the same volume of material, compared to a pure edge dislocation. Therefore, the trapped dislocation dipole is broken at the corners where the climb rate is higher than at the middle, turning the edge dipole into an elongated interstitial prismatic loop and two helices, as shown in gure 7(d). As one may expect, the elongated loop would continue to evolve into small loops by sequential pinching off as discussed in subsection 3.2.1. This allows further glide along the direction parallel to the Burgers vector. This mechanism of dipole breakup is believed to play a signi cant role in releasing the strain-hardening in the subsequent annealing of deformed crystals [13]. The remaining helices may shrink or expand perpendicular to their Burgers vector by interacting with prismatic loops or by bulk-diffusion controlled climb [38,77].

Conclusions
Dislocations can provide short circuit diffusion paths for atoms. This fast atomic transport, along the dislocation core region, known as core diffusion, accelerates the diffusion by more than three orders of magnitude compared to lattice diffusion. This allows dislocation motions perpendicular to the original slip system, known as dislocation self-climb (conservative climb) and is of particular importance in low-temperature creep and post-irradiation annealing.
A variational principle is presented for the analysis of problems in which fast dislocation core diffusion is the dominant mechanism for material redistribution. A new nite element discretisation method is employed to accelerate the simulation of dislocation self-climb. We have demonstrated here, that the paired-linear element based method has a much higher computational ef ciency compared with the previous linear element based method [1]. This becomes more evident as the segment number increases. The acceleration in computation enables the simulation of collective dislocation motion, including both glide and self-climb, based on the current method.
An analytical solution for the evolution of an isolated elliptical prismatic loop was derived based on the variational principle, and compared with the numerical method. The analytical solution provides a good approximation for predicting the evolution of elliptical loops with small aspect ratios. However, as the aspect ratio increases, the high curvature at the two end of the ellipse cause the loop to evolve into a bone-like shape initially, and gradually converge into a circular loop, which the analytic solution fails to predict. This interesting discrepancy leads to the study of edge dipoles or elongated prismatic loops, for which the aspect ratios are much higher. The evolution of two representative edge dipoles, the close-ended jog dragging dipole and the open-ended edge trapping dipole, were examined with the new DDD method. Results show that a string of prismatic loops is formed by the sequential pinchingoff of the dipole, rather than spontaneous break-up, which agrees well with experimental observations.