A phase field crystal theory of the kinematics of dislocation lines

We introduce a dislocation density tensor and derive its kinematic evolution law from a phase field description of crystal deformations in three dimensions. The phase field crystal (PFC) model is used to define the lattice distortion, including topological singularities, and the associated configurational stresses. We derive an exact expression for the velocity of dislocation line determined by the phase field evolution, and show that dislocation motion in the PFC is driven by a Peach-Koehler force. As is well known from earlier PFC model studies, the configurational stress is not divergence free for a general field configuration. Therefore, we also present a method (PFCMEq) to constrain the diffusive dynamics to mechanical equilibrium by adding an independent and integrable distortion so that the total resulting stress is divergence free. In the PFCMEq model, the far-field stress agrees very well with the predictions from continuum elasticity, while the near-field stress around the dislocation core is regularized by the smooth nature of the phase-field. We apply this framework to study the rate of shrinkage of an dislocation loop seeded in its glide plane.


Introduction
Plasticity in crystalline solids primarily refers to permanent deformations resulting from the nucleation, motion, and interaction of extended dislocations. Classical plasticity theories deal with the yielding of materials within continuum solid mechanics [Hill, 1998, Wu, 2004. Deviations from elastic response are described with additional variables (e.g., the plastic strain), which effectively describe the onset of plasticity (yielding criteria), as well as the mechanical properties of plastically deformed media (e.g., work hardening). A macroscopic description of the collective behavior of dislocation ensembles is thus achieved, usually assuming homogeneous media for large systems. In crystal plasticity, inhomogeneities and anisotropies are accounted for, with the theory having been implemented as a computationally efficient finite element model [Roters et al., 2010, Pokharel et al., 2014. These theories are largely phenomenological in nature, and rely on constitutive laws and material parameters to be determined by other methods, or extracted from experiments. They can be finely tuned, but sometimes fail in describing mesoscale effects [Rollett et al., 2015]. On the other hand, remarkable mesoscale descriptions have been developed by tracking single dislocations [Kubin et al., 1992, Bulatov et al., 1998, Sills et al., 2016, Koslowski et al., 2002, Rodney et al., 2003. These approaches typically evolve dislocation lines through Peach-Koehler type forces while incorporating their slip system, mobilities, and dislocation reactions phenomenologically. Stress fields are described within classical elasticity theory [Anderson et al., 2017]. Since linear elasticity predicts a singular elastic field at the dislocation core, theories featuring its regularization are usually exploited. Prominent examples are the non-singular theory obtained by spreading the Burgers vector isotropically about dislocation lines [Cai et al., 2006], and the stress field regularization obtained within a strain gradient elasticity framework [Lazar and Maugin, 2005]. Plastic behavior then emerges when considering systems with many dislocations and proper statistical sampling [Devincre et al., 2008]. Still, the accuracy and predictive power of these approaches depend on how well dislocations are modeled as isolated objects. In this context, mesoscale theories that require a limited set of phenomenological inputs are instrumental in connecting macroscopic plastic behavior to microscopic features of crystalline materials.
The Phase Field Crystal (PFC) model is an alternative framework to describe the nonequilibrium evolution of defected materials at the mesoscale [Elder et al., 2002, Emmerich et al., 2012, Momeni et al., 2018. Within the phase q (4) q (5) q (6) q (7) q (8) q (9) q (10) q (11) q (12) q x q y q z Figure 1: (a) A dislocation line C consisting of points r characterized by the tangent vector t and the Burgers vector b at that point. The difference r − r from a point r to a point r on the line can be decomposed into a 2D in-plane vector ∆r ⊥ , which is the projection of r − r onto the plane N normal to t and a distance |∆r | from this plane. In this figure, ∆r · t = −3.47a 0 and |∆r ⊥ | = 3.42a 0 . (b) The N = 12 primary reciprocal lattice vectors {q (n) } 12 n=1 of length q 0 of a bcc lattice (Eq. (20)). Higher modes (dots) correspond to higher harmonics {p n } n>N in the expansion of the equilibrium phase-field configuration ψ eq , Eq. (3). lines in three dimensions down to the defect core. In the case of a point defect, it was shown in Ref. [Salvalaglio et al., 2020] that the stress field at the core agrees with the predictions of the non-singular theory of Ref. [Cai et al., 2006], and with gradient elasticity models Maugin, 2005, Lazar, 2017], indicating that the results obtained here can serve as benchmarks for similar theories in three dimensions. The specific example of a dislocation loop in a bcc lattice is considered, and the far-field stresses given by the ψ field are shown to coincide with predictions by continuum elasticity.
The rest of the paper is structured as follows. In Sec. 2, we introduce the theoretical method used to define topological defects from a periodic ψ-field. This allows us to define a dislocation density tensor in terms of the phase field (Eq. (12)), and obtain the dislocation line velocity (Eq. (16)). These are key results, which are applied in several examples in Sec. 3. First, we use the PFC model to numerically study the shrinkage of a dislocation loop in a bcc lattice. Then, we show analytically that Eq. (16) captures the motion of dislocations driven by a Peach-Koehler type force, and hence by a local stress. Finally, we introduce the PFC-MEq model, and compare the shrinkage of the dislocation loop under PFC and PFC-MEq dynamics. While the results are qualitatively similar for the case of a shear dislocation loop, the constraint of mechanical equilibrium causes the shrinkage to happen much faster. We finally confirm that the stress field derived from the ψ field in the PFC-MEq model agrees with that which would follow from continuum elasticity theory, with the same singular dislocation density as source, and with no adjustable parameters.

Kinematics of a dislocation line in three dimensions
Dislocations in 3D crystals are line defects, where each point r on the line C is characterized by the tangent vector t at that point and a Burgers vector b, see Figure 1(a). By introducing a local Cartesian plane N normal to t , the distance of an arbitrary point r to a point r on C can be decomposed into an in-plane vector ∆r ⊥ ⊥ t and a vector ∆r t , i.e. r − r = ∆r ⊥ + ∆r . A deformed state can be described by a displacement field u and, in the presence of a dislocation, u is discontinuous across a surface (branch cut) spanned by the dislocation, given by where u + and u − are the values of the displacement field at each side of the branch cut, respectively. We use the negative sign convention relating the contour integral with the Burgers vector. Here, Γ is a small circuit enclosing the dislocation line in the N -plane, directed according to the right-hand rule with respect to t . The dislocation density tensor associated with the line is [Lazar, 2014] where b j is the j component of the Burgers vector of the line, and dl i = t i dl is the line element in the direction of the line. δ (2) i (C) is a short-hand notation for the delta function, with dimension of inverse area, locating the position of the dislocation line for each component i of the dislocation density tensor. It is defined by the line integral over the dislocation line of the full delta function (which scales as inverse volume). The dislocation density tensor is defined so that N d 2 r ⊥ α i j t i = b j , where we are using the Einstein summation convention over repeated indices.
In the PFC models, a crystal state is represented by a periodic phase field ψ(r) of a given crystal symmetry. A reference crystalline lattice, is defined by a set of N primary (smallest) reciprocal lattice vectors {q (n) } N n=1 of length q 0 , and higher harmonics {p n } n>N , also on the reciprocal lattice but with |p n | > q 0 (see, e.g., {q (n) } with |q (n) | = q 0 for a bcc lattice in Fig. 1(b)). The lattice constant of the crystal is then given by a 0 ∼ 2π/q 0 . This represents a perfect crystal configuration in the absence of defects and distortion, where the average valueψ and the amplitudes η n are constants. In the phase-field crystal theory presented in Refs. [Elder et al., 2002, Elder andGrant, 2004], near the solid-liquid transition point, only the terms from the primary reciprocal lattice vectors contribute to ψ eq , while in general for more sharply peaked density profiles, there are also contributions from the higher order harmonics {p n } n>N . For a distorted crystal lattice, the mode amplitudes η n become complex scalar fields, henceforth named complex amplitudes η n (r), such that In this section, we provide an accurate description of dislocation lines as topological defects in the phase of the complex amplitudes η n (r). We generalize the method of tracking topological defects as zeros of a complex order parameter as introduced in Refs. [Halperin, 1981] and [Mazenko, 1997], and apply it to accurately derive the kinematics of dislocation lines. Given a phase field configuration ψ(r), the complex amplitudes can be found by a demodulation as described in Appendix A.1. Decomposing each amplitude η n (r) = ρ n (r)e iθ n (r) , into its modulus ρ n (r) and phase θ n (r), we have that for a perfect lattice, θ (0) n = 0 and ρ n is constant. Displacing a lattice plane by a slowly varying u transforms the phase as θ n → θ (0) n − q (n) · u. Thus, the phase provides a direct measure of the displacement field u(r) relative to the reference lattice, i.e., where q (n) i denotes the i-th Cartesian coordinate of q (n) . It is possible to invert Eq. (5), and solve for the displacement field u as function of the phases θ n and reciprocal vectors. We use the following identity which is valid for lattices with cubic symmetry, where all primary reciprocal lattice vectors have the same length q 0 (see Appendix B) so that the displacement u is given by Eq. (7) shows that a dislocation line, which introduces a discontinuity in the displacement field, leads to a discontinuity in the phases θ n (r). This is the first key insight, which we illustrate in Fig. 2. By using Eq. (5) and the fact that the Burgers vector b is constant along the dislocation line, we relate the Burgers vector to the phase θ n as where s n is the (integer) winding number of the phase θ n around the dislocation line. That s n is an integer follows from the fact that while θ n (r) may have a discontinuity across the branch cut, the complex amplitude η n (r) is well-defined and continuous everywhere. Therefore the circulation of the phase must be an integer multiple of 2π. By the same reasoning, for an amplitude for which q (n) · b 0, at the dislocation line, the phase θ n (r) is undefined (singular), so the modulus ρ n (r) must go to zero for η n (r) to remain continuous. This is the second key insight, which allows us to identify the location of the dislocation line with the zeros of the complex amplitudes η n (r). The complex amplitude η n (r) is isomorphic to a 2-component vector field Ψ(r) ≡ (Ψ 1 (r), Ψ 2 (r)) = ( (η n (r)), (η n (r))). The study of how to track zeros of any dimensional vector field in any dimensions was introduced in Ref. [Halperin, 1981]. The orientation field Ψ(r)/|Ψ(r)| is continuous wherever |Ψ(r)| 0 and supports 1D topological defects in 3 dimensions which are located precisely where |Ψ(r)| = 0. The topological line density ρ i of the line C, which satisfies d 2 r ⊥ ρ i = s n t i , is given by ρ = s n δ (2) (C) ρ i = s n δ (2) i (C) .
Like δ (2) i (C), the dimension of ρ i is that of a two-dimensional vector density. This topological charge density is expressed explicitly in terms of the real-valued positions C = {r } of the topological defect line. Since these positions coincide with the zero-line of the vector field Ψ(r), it is possible to relate the expression to the delta-function locating the zeros of Ψ(r), through the transformation law s n δ (2) (C) = D i (r)δ (2) (Ψ(r)), with the determinant vector field D i (r) = i jk (∂ j Ψ 1 (r))(∂ k Ψ 2 (r)). Comparing this to Eq. (2), using Eq. (8) and re-expressing D (n) i (r) (with the added superscript n) in terms of the complex amplitude η n (r), we end up with the central equation for tracking the evolution of the dislocation density where δ (2) (η n ) = δ( (η n ))δ( (η n )) and D (n) (r) = ∇ (η n (r)) × ∇ (η n (r)) D (n) i (r) = i jk (∂ j (η n (r)))(∂ k (η n (r))) .
In the following, for ease of notation, we suppress the explicit positional dependence of α i j , D (n) i and η n . The dislocation line is located at η n = 0, which is the intersection of the surfaces (η n ) = 0 and (η n ) = 0. As we see from its definition, D (n) is perpendicular to both these surfaces and is thus directed along the tangent to the line. We can reconstruct the dislocation density tensor from an appropriate summation over the modes with singular phases, namely by multiplying Eq. (10) by q (n) j , summing over the reciprocal modes and using Eq. (6) to arrive at Having a closed form of the dislocation density in terms of the complex amplitudes η n , we now turn to deriving a closed form expression for its kinematic in terms of the time evolution of η n . Taking the time derivative of Eq.
(2), we show in Appendix C.1 that for a dislocation density tensor described by a single loop or string, we have and V is a vector field defined on the string by the velocity of the line segment perpendicular to the tangent vector.
Taking the time derivative of Eq. (12), we show in Appendix C.2 that we get and . Note that J l j depends on ∂ t η n , and hence on the law governing the temporal evolution of the phase field. J (α) l j is the well-known expression in terms of the dislocation velocity and J l j is what we predict from the evolution of the phase field crystal density ψ. Under the assumption that both currents are equal, we show in the following that we are able to determine the dislocation velocity directly from the evolution of the phase field ψ at the dislocation core. We have checked numerically that the dislocation velocity predicted with this assumption is in excellent agreement with the one computed by tracking the position of the dislocation line at successive time steps.
By contracting Eq. (10) with D (n) i , we can express the delta-function in terms of the dislocation density tensor , which we can insert into Eq. (14). Then, by equating J l j and J (α) l j at a point r on the dislocation line, where α ik = t i b k δ (2) (∆r ⊥ ), we get after contracting with b and integrating the delta-functions in N (details in Appendix D) where v is the velocity of the dislocation node at r . Since t ⊥ v , we can easily invert this relation to find v , and using that Eqs. (12) and (16) are the key results of this paper. Eq. (12) defines the dislocation density tensor from the demodulated amplitudes η n of the phase field, while Eq. (16) gives an explicit expression for the dislocation line velocity. Both equations bridge the continuum description of the dislocation density and velocity with the microscopic scale of the phase field.

Dislocation motion in a bcc lattice
We apply here the framework developed in Sec. 2 to a phase field crystal model of dislocation motion in a bcc lattice [Elder et al., 2002, Elder and Grant, 2004, Emmerich et al., 2012. The free energy F ψ is a functional of the phase field ψ over the domain Ω, given by where L = q 2 0 + ∇ 2 , and ∆B 0 , B x 0 , V, and T are constant parameters . The dissipative relaxation of ψ reads as ∂ψ ∂t = Γ∇ 2 δF ψ δψ .
with constant mobility Γ. We will refer to Eq. (18) as the "classical" PFC dynamics. As a characteristic unit of time given these model parameters, we use τ = (ΓB x 0 q 6 0 ) −1 . For appropriate parameter values, the ground state of this energy is a bcc lattice which is well described in the one mode approximation where ψ 0 is the average density, η 0 is the equilibrium amplitude found by minimizing the free energy (Eq. (17)) with this ansatz for ψ(r), and {q n } are the N = 12 smallest reciprocal lattice vectors with q (n) = −q (n−6) for n = 7, .., 12, see Fig. 1(b). Figure 3 shows one bcc unit cell of a phase-field initialized in the one-mode approximation. Given the equilibrium configuration, the lattice constant a 0 will be used as the characteristic unit of length and the shear modulus µ calculated from the phase-field will serve as the characteristic unit of stress [Skogvoll et al., 2021a]. As we see, the functional form of the free energy determines the base vectors q (n) , and no further assumptions about slip systems or constitutive laws for dislocation velocity (or plastic strain rates) need to be introduced.
The model parameters (∆B 0 , B x 0 , T, V, and Γ) and variables (F ψ , ψ, r, and t) can be rescaled to a dimensionless form in which B x 0 = V = q 0 = Γ = 1, thus leaving only three tunable model parameters: the quenching depth ∆B 0 , T and the average density ψ 0 (due to the conserved nature of Eq. (18)). All simulations are performed in these dimensionless units as described in Sec. Appendix A.3.

Numerical analysis: shrinkage of a dislocation loop
In order to have a lattice containing one dislocation loop as the initial condition, we consider first the demodulation of the ψ field in the one mode approximation. A dislocation loop is introduced into the perfect lattice by multiplying the equilibrium amplitudes by complex phases η 0 → η n (r) with the appropriate charges s n (see Appendix A.4) and then reconstructing the phase field ψ through Eq. (19). We then integrate Eq. (18) forward in time as detailed in Appendix A.3. A fast relaxation follows from the initial configuration with the loop. This relaxation leads to the regularization of the singularity at the dislocation line (η n → 0 for s n 0) as achieved in PFC approaches [Skaugen et al., 2018a, Salvalaglio et al., 2019. From then onward, ψ evolves in time leading to the motion of the dislocation line which may be analyzed by the methods outlined in Sec. 2, using the amplitudes {η n } extracted from ψ extracted as detailed in Appendix A.1.
Numerically, we approximate the delta function in Eq. (12) as a sharply peaked 2D Gaussian distribution, i.e., δ (2) (η n ) exp(− |η n | 2 2ω 2 )/(2πω 2 ) with a standard deviation of ω = η 0 /10. Near the dislocation line, the dislocation density α i j thus takes the form of a sharply peaked function, which can be treated numerically. The decomposition of α i j into its outer product factors t i and a Burgers vector density B j = b j δ (2) (∆r ⊥ ) is done by singular value decomposition (see Sec. Appendix A.2), and the Burgers vector of the point is extracted by performing a local surface integral in N . We prepare a 35 × 35 × 35 unit cell 3D PFC lattice on periodic boundary conditions with a resolution of ∆x = ∆y = ∆z = a 0 /7. A dislocation loop is introduced as the initial condition in the slip system given by a plane normal [−1, 0, 1] with slip direction (Burgers vector) a 0 2 [1, −1, 1]. Figure 4(a) shows the initial dislocation density decomposed as described, where we also have calculated the velocity v at each point given by Eq (16).
In order to obtain the velocity of the dislocation loop segments, we identify M nodes on the loop and evaluate Eq. (16) by using numerical differentiation of the ψ field to calculate the amplitude currents J (n) l . To serve as a benchmark, we also calculate the circumference l C of the dislocation loop C at each time (further details in Appendix A.5), so that we compare the rate of shrinkage |∂ t R| (solid blue line in Fig. 4 (dashed blue line in Fig. 4(b)) where v (m) is the velocity of the dislocation line at node m, calculated by the velocity formula Eq. (16). |∂ t R| andv should agree in the case of the shrinking of a perfectly circular loop and the figure shows excellent agreement between the two. Interestingly, we observe that both are sensitive to the Peierls like barriers during their motion, as shown by the oscillations in Fig. 4(b)). The maxima are separated by 2πa 0 , confirming that the oscillation is related to the motion of a loop segment over one lattice spacing a 0 [Boyer and Viñals, 2002]. This observation confirms that even though Eqs. (12) and (16) are continuum level descriptions of the system, they still exhibit behavior related to the underlying lattice configuration. The initial fast drop in velocity is due to the fast relaxation of the initial condition. The evolution of the variables under the dynamics of Eq. (18) are shown together with the evolution given by the PFC-MEq model which will be introduced in Sec. 3.3.

Theoretical analysis: Peach Koehler law
In this section, we show that the general expression Eq. (16) of the defect velocity agrees with the dissipative motion of a dislocation as given by the classical Peach-Koehler force [Pismen, 1999, Kosevich, 1979. To calculate an analytical expression for the amplitude currents J l , we employ the amplitude formulation of the PFC model, which directly expresses the free energy and dynamical equations in terms of the complex amplitudes η n [Goldenfeld et al.,  2005, Athreya et al., 2006, Salvalaglio andElder, 2022]. For our lattice symmetry, real valuedness of ψ requires that η n+6 = η * n , and the dynamical equations need only consider the amplitudes {η n } 6 n=1 . By substituting Eq. (19) in F ψ and integrating over the unit cell, under the assumption of slowly-varying amplitudes, one obtains the following free energy as a function of the complex amplitudes, where G n = ∇ 2 +2iq n ·∇ and Φ = 2 6 n=1 |η n | 2 . f s ({η n }, {η * n }) is a polynomial in η n and η * n that depends in general on the specific crystalline symmetry under consideration [Goldenfeld et al., 2005, Salvalaglio and Elder, 2022 (here bcc, see Appendix E for its expression). Equation (23) is obtained when considering a set of vectors q of length q 0 , while similar forms may be achieved when considering different length scales , Salvalaglio et al., 2021. The evolution of η n , which follows from Eq. (18) is [Goldenfeld et al., 2005, Salvalaglio andElder, 2022], with where the last term comes from the nonlinear contributions ψ 3 and ψ 4 in the local free energy density, and depend on the other amplitudes {η m } m n . However, for the amplitudes that go to zero at the defect, it can be shown that ∂ f s ∂η * n = 0 at the defect (for more details, see Appendix E). Thus, the evolution of η n near the defect core is dictated solely by the non-local gradient term, namely Furthermore, this implies that the complex amplitude η n of a stationary defect satisfy G 2 n η (0) n = 0 at the core. We now add an imposed, smooth displacementũ to the amplitudes as η n = η (0) n e −iq n ·ũ to represent the far-field displacement induced by a different line segment, defect, or externally applied loads [Skaugen et al., 2018a]. This displacement is in addition to the discontinuous displacement field u, described in Sec. 2, which is captured by stationary solution η (0) n and defines the Burgers vector of the dislocation line (Fig. 2). Inserting this ansatz of the complex amplitudes into Eq. (11), and in the approximation of small distortions, |∇ũ| 1, we find where D (n),0 i is the determinant vector field calculated from η (0) n . The corresponding defect density current is Arguably, the simplest solution of Eq. (26) is the isotropic, simple vortex η (0 ) n which is linear with the distance from the core and s n = ±1. At a node r on the dislocation line, η (0 ) can be written in terms of the Cartesian coordinates x ⊥ , y ⊥ in the plane N (Sec. 2), where it takes the form η (0 ) n = κ(x ⊥ + is n y ⊥ ), with κ a proportionality constant. The gradients of η (0 ) n can be evaluated in these coordinates and gives at r , , from which we get the current in terms of the local tangent vector t . At r , we also get D (n) i = κ 2 s n t i , which leads to an expression of the dislocation velocity (where the proportionality constant κ cancels out), given by whereσ l j is the stress tensor for a bcc PFC that has been deformed byũ [Skogvoll et al., 2021a], Thus, the velocity of the dislocation line is proportional to the stress on the line. In vectorial form, this equation reads with isotropic mobility M = Γπ/(|b| 2 η 2 0 ). A stationary dislocation induces a stress field σ (0) i j , but only the imposed stressσ i j appears in the equation above. This is analogous to how the stress field of the dislocation itself is not included when the Peach-Koehler force as calculated [Kosevich, 1979]. Thus, if σ ψ i j is the configurational stress of the phase field at any given time, the part responsible for dislocation motion is the imposed stress Note that the stationary solution necessarily satisfies mechanical equilibrium, ∂ j σ (0) i j , so that if the configurational PFC stress σ ψ i j is in mechanical equilibrium, so is the imposed stressσ i j on the dislocation segment. The imposed stress used can be attributed to external load, other dislocations, or other parts of the dislocation loop. The framework predicts a defect mobility which is isotropic and does not discriminate between dislocation climb and glide motion. Numerically however, we have seen that at deeper quenches ∆B 0 , climb motion is prohibited in the PFC model. The result in this section should therefore be interpreted as a first-order approximation, valid at shallow quenches. This apparent equal mobility for glide and climb may result from the employment of the amplitude phase-field model (which is only exact for |∆B 0 | → 0) or the assumption of an isotropic defect core in the calculation.

PFC dynamics constrained to mechanical equilibrium (PFC-MEq)
In the previous section, we found that the motion of a dislocation is governed by a configurational stress σ ψ i j which derives from the PFC free energy. Since this stress is a functional only of the phase field configuration, it does not satisfy, in general, the condition of mechanical equilibrium. References [Skaugen et al., 2018a, Skogvoll et al., 2021a give an explicit expression for this stress defined as the variation of the free energy with respect to distortion, where · is a spatial average over 1/q 0 in order to eliminate the base periodicity of the phase field (see Appendix A.1).
In this section, we discuss a modification of the PFC in three dimensions and in an anisotropic lattice so as to maintain elastic equilibrium in the medium while ψ evolves according to Eq. (18). Let ψ (U) be the field that results from the evolution defined by Eq. (18) alone. At each time, we define where u δ is a small continuous displacement computed so that the configurational stress associated with ψ(r) is divergence free. We now show a method to determine u δ . Suppose that at some time t the PFC configuration ψ has an associated configurational stress σ ψ,U i j (from Eq. (34), where ∂ j σ ψ,U i j 0). Within linear elasticity, the stress σ ψ i j after displacement of the current configuration by u δ is given by where C i jkl is the elastic constant tensor, and e δ i j = 1 2 (∂ i u δ j + ∂ j u δ i ). u δ is determined by requiring that By using the symmetry i ↔ j of the elastic constant tensor, we can rewrite this equation explicitly in terms of u δ , where is the body force from the stress [Skogvoll et al., 2021a]. The quantity f is the free energy density from Eq. (17). Given the periodic boundary conditions used, the system of equations (38) is solved by using a Fourier decomposition with the Green's function for elastic displacement in cubic anisotropic materials [Dederichs and Leibfried, 1969]. Once u δ is obtained, ψ is updated according to Eq. (35), and evolved according to Eq. (18) from its current state ψ(t) to ψ (U) (t + ∆t). Note that Eqs. (38) can, in general, be solved for any elastic constant tensor, so that the method introduced is not limited to cubic anisotropy. Since the state ψ (U) can only be updated according to Eq. (35) every ∆t, this effectively sets a time scale of elastic relaxation in the model. We found that the numerical discretization scheme for imposing mechanical equilibrium at every ∆t has a slow convergence with decreasing time resolution. Thus, the rate of loop shrinkage also depends slightly on ∆t. This is further discussed in Appendix A. Figure 4 contrasts numerical results for the evolution of an initial dislocation loop with and without using the method just described. The computed line velocities are very different as they are highly sensitive to the local stress experienced by the dislocation loop segments. This stems from the fact that under classical PFC dynamics, the stress is always given by σ ψ,U i j , and a consequence of the results from Sec. 3.2 is that the velocity of an element of the defect line will be quite different depending on whether the stress acting on it is σ ψ,U i j or σ ψ i j . Figure 5 shows the dislocation loop after its circumference has shrunk to 90% of its initial value, and the resulting xz component of the stress for both models. As expected, the correction provided by the PFC-MEq model is necessary to relax the stress (a) (b) Figure 5: In-plane sections (y = 17.5a 0 ) of the configurational stress σ ψ xz /µ for the dislocation loop after shrinking to 90% of its initial circumference under (a) PFC dynamics and (b) PFC-MEq dynamics. Because the latter evolves faster, the snapshots are taken at different times, namely t = 389.0τ and t = 34.4τ, respectively. A lot of residual (unrelaxed) stress is visible in the configurational stress for the classical PFC model. originating from the initial loop. The figure shows a large residual stress far from the dislocation loop that can only decay diffusively in the standard phase field model. Indeed, we have verified numerically that the configurational stress is only divergence-less for the PFC-MEq model. We note that in our set the loop is seeded in a glide plane, thus its shape remains approximately circular for both models, while the shrinkage rate is different. Note that with the addition of this advection step, the model is no longer guaranteed to be fully dissipative.
The problem addressed in this section involves finding the elastic distortion u kl (which away from defects it can be written as u kl = ∂ k u l for a displacement field u) given the dislocation density tensor α i j as a state variable [Acharya et al., 2019]. The first part is the incompatibility of the elastic distortion and the second is the mechanical equilibrium condition on u kl where C i jkl is the tensor of elastic constants, and (S ) denotes the symmetric part of the tensor. Equation (40) has a non trivial kernel consisting of gradients of vector fields ∇u δ . This vector field is determined by Eq. (41) given appropriate boundary conditions that guarantee uniqueness. A computational method for solving for u kl and u δ , using the dislocation density as a state variable, was first given in Ref. [Roy and Acharya, 2005]. The main difference between this reference and the method outlined in this section is that, since the incompatibility of the distortion is captured by the state of the phase field, we only need to solve for the compatible part of the distortion using the force density g ψ from the phase field as a source.
While the stress profile shown in Fig. 5(b), can be shown numerically to have vanishing divergence, we would like to see a direct comparison of the stress with the prediction from continuum elasticity. As the model purports to evolve the phase-field at mechanical equilibrium, and we are able to extract the dislocation density from the phase-field at any time through Eq. (12), this amounts to the problem of finding the stress tensor for a given dislocation density, under the constraint of mechanical equilibrium and with periodic boundary conditions (zero surface traction). This problem was adressed in Ref. [Brenner et al., 2014], and in Appendix A.6, we show how we solve Equations (40-41) to derive the equilibrium stress field from α i j using spectral methods. Figure 6 shows all the stress components after the dislocation loop has shrunk to 90% of its initial diameter for both dynamical models, as well as the stress σ (α) computed directly from the dislocation density tensor. 1 Note that the mean value of the components of σ (α) i j , is not determined by Eqs. (40) -(41), and is set to zero. In this comparison, we have also subtracted from σ ψ i j its mean value. As expected, the stresses obtained from the PFC-MEq model agree well with σ (α) i j . The small differences observed are due to the fact that the configurational stress determined by ψ is naturally regularized by the lattice spacing and the finite defect core, whereas the stress σ (α) i j is for a continuum elastic medium with a singular dislocation source (numerically, the δ-functions in Eq. (12) is regularized by an arbitrary width of the Gaussian approximation). Investigating exactly which length scale of core regularization derives from the PFC model is an open and interesting question that we will address in the future.

Conclusions
We have introduced a theoretical method, and the associated numerical implementation, to study topological defect motion in a three dimensional, anisotropic, crystalline PFC lattice. The dislocation density tensor and velocity are directly defined by the spatially periodic phase field, where dislocations are identified with the zeros of its complex amplitudes.
To illustrate the method, we have studied the motion of a shear dislocation loop, and found that it accurately tracks the loop position, circumference, and velocity. As an application, we have shown that under certain simplifying assumptions, the overdamped dislocation velocity follows from the Peach-Koehler force, with the defect mobility determined by equilibrium lattice properties. We have introduced the PFC-MEq model for three dimensional anisotropic media which constrains the classical PFC model evolution to remain in mechanical equilibrium, and shown that loop motion is much faster with this modification. The PFC-MEq model produces stress profiles that are in agreement, especially far from the defect core, to stress fields directly computed from the instantaneous dislocation density tensor.
In summary, we have presented a comprehensive framework, based on the phase field crystal model for the analysis of dislocation motion in crystalline phases in three spatial dimensions. Starting from a free energy that has a ground state of the proper symmetry, the model naturally incorporates defects, the associated topological densities, and the resulting defect line kinematic laws that are compatible with topological density conservation. Configurational stresses induced by defects are defined and analyzed, and shown to lead to a Peach-Koehler type force on defects, with an explicit expression for the line segment mobility given. we can find the amplitudes using the principle of resonance under coarse graining. Coarse grainingX with respect to a length scale a 0 is introduced as a convolution with a Gaussian filter function Given the PFC configuration of Eq. (A.1), to find η n (r), we multiply by e −q (n) ·r and coarse grain to get ψ(r)e −q (n) ·r =ψ(r) e −q (n) ·r + n η n (r) e i(q (n ) −q (n) )·r = η n (r), where we have used the slowly varying nature of the complex amplitudes to pull them out of the coarse graining operation and used the resonance condition e i(q (n ) −q (n) )·r = δ nn [Skogvoll et al., 2021a].
Appendix A.2. Dislocation density tensor decomposition A singular value decomposition of α is introduced as α = UΣV T , where Σ is a diagonal matrix containing the singular values of α, and U and V are unitary matrices containing the normalized eigenvectors of (αα T ) and (α T α), respectively. We assume that the dislocation density tensor can be written as the outer product of the unitary tangent vector t and a local spatial Burgers vector density B(r), i.e., α i j = t i B j . Under this assumption, one finds Σ with only one non zero singular value, |B|, and the columns of U and V that correspond to this singular value will be t and B/|B|, respectively.

Appendix A.3. Evolution of the phase field
The dimensionless parameters for the bcc ground state are set to: ∆B 0 = −0.3, T = 0 and ψ 0 = −0.325. Lengths have been made dimensionless by choosing |q (n) | = q 0 = 1, yielding a bcc lattice constant a 0 = 2π √ 2. In all simulations, the computational domain is given by 35 × 35 × 35 base periods of the undistorted bcc lattice, with grid spacing ∆x = ∆y = ∆z = a 0 /7. Periodic boundary conditions are used throughout. Equation (18) is integrated forward in time with an explicit method [Cox and Matthews, 2002], and ∆t = 0.1. A Fourier decomposition of the spatial fields is introduced to compute the spatial derivatives of the fields, while nonlinear terms are computed in real space. N is the plane normal to the tangent vector t upon which we impose a Cartesian coordinate system to determine the angles θ 1 , θ 2 that are used to construct the (inset) initial amplitude phase configuration. For more details, see Section Appendix A.4.

Appendix A.3.1. Mechanical equilibrium
We implement the correction scheme of Eq. (35) between every timestep ∆t. If u max = max r∈Domain (u δ (r)) > 0.1a 0 , we rescale u δ so that u max = 0.1a 0 , and repeat the process again until elastic equilibrium is achieved. Typically, when initializing the PFC field with a dislocation, around 5 such iterations are needed, after which, u max is on the order of 0.01a 0 at each correction step.
The dislocation loop shrink velocity is sensitive to the time interval ∆t between each equilibration correction. As shown in Fig. 4, the effect of imposing this correction at every time interval ∆t = 0.1 accelerates the annihilation process by approximately a factor of |v PFCMEq,∆t=0.1 |/|v PFC | ≈ 7.5. A slow convergence in the limit ∆t → 0 is observed, where we have estimated that the shrink velocity increases up to |v PFCMEq,∆t→0 |/|v PFC | ≈ 9.8. However, to reach this numerical convergence is computationally demanding. Indeed, this slow convergence suggests that the time scale of the elastic field relaxation is important for the process of shear dislocation loop shrinkage. For static problems however, such as obtaining regularized stress profiles for dislocation loops, or defect nucleation under quasi-static loading, this slow convergence is not an issue.

Appendix A.4. Initializing a dislocation loop in the PFC model
In this section, we show how to multiply the initial amplitudes η 0 with complex phases, to produce a dislocation loop with Burgers vector b in a slip plane given by normal vector n (see Sec. 3.1). Given a point r, it belongs to a plane N perpendicular to t for some point r on the dislocation loop (see Figure A.7). This plane also intersects the diametrically opposed point r of the dislocation loop. If r 0 is the center of the loop, the distance vector r − r 0 lies in N . Let (m 1 , m 2 ) be the first and second coordinate in the Cartesian coordinate system defined by the right-handed orthonormal system {(n × t ), n, t } centered at r 0 . If m 1 > 0, we get from geometrical considerations Both m 1 and m 2 are thus determined by r, r 0 and the normal vector to the loop plane n. θ 1 (θ 2 ) is the angle between r − r (r − r ) and n × t in the plane N and are found numerically by using the four-quadrant inverse tangent atan2(y, x), so that θ 1 = atan2 (m 2 , m 1 + R) (A.6) where R is the radius of the loop. For each point r, we determine θ 1 (r) and θ 2 (r) according to the equations above and initiate the PFC with the phases η n = η 0 e is n (θ 1 (r)−θ 2 (r)) , (A.8) where s n = 1 2π q (n) · b is given in Table A To calculate numerically the perimeter of a dislocation loop, recall that where we have added a subscript (l) onto r to emphasize that it is the point on the loop as indexed by the line element dl. Taking the double dot product with itself, we find The contributions to this integral will only come from points on the loop C and only when r (l) = r (m) , where dl i = dm i , so dl i dm i = (dl) 2 = |dl i | 2 = |dl i ||dm i |. Thus Taking the square root and integrating over all space, we find where L is the perimeter of the dislocation loop. Thus, Appendix A.6. Direct computation of stress fields The dislocation density tensor is calculated directly from the phase field ψ through Eq. (12). The general method of solving Eqs. (40-41) on a periodic medium is given in Ref. [Brenner et al., 2014] given α i j , where also the uniqueness of the elastic fields is proven given appropriate conditions on the dislocation density α i j . In the present case, the conditions on α i j are automatically satisfied as it is calculated from the phase-field. In this section, we thus show for our computational setup, how we compute the Green's function in the relating the distortion u i j to the dislocation density tensor α i j as a source. Since (40-41) given the periodic boundary conditions can be solved uniquely, we Fourier transform both sets of equations and add the condition of mechanical equilibrium (Eq. (41)) to the diagonal equations (i = k) in Eq. (40), which gives in Fourier space where there is no summation over (i), and we have multiplied the elastic constant tensor by i/µ where µ is the shear modulus of the cubic lattice, and C i jkl = λδ i j δ kl + µ(δ ik δ jl + δ il δ jk ) + γδ i jkl . By defining the 1D vectorsŨ andα as U T = u 11 , u 12 , u 13 , u 21 , u 22 , u 23 , u 31 , u 32 , u 33 , α T = α 11 , α 12 , α 13 , α 21 , α 22 , α 23 , α 31 , α 32 , α 33 , we rewrite Eq. (A.14) more compactly as where the explicit form of M(q) in the case of cubic anisotropy is given by M(q) can be inverted to yield the Fourier transform of the distortion U, OnceŨ (denoted byũ kl in components) is known, we compute the stress field in mechanical equilibrium The dislocation density α ik as obtained from the phase field as in Eq. (12) has a very small divergence due to numerical round-off errors. We impose ∂ i α ik = 0 explicitly before evaluatingσ, which improves numerical stability.

Appendix B. Inversion formula for highly symmetric lattice vector sets
In inverting Eq. (5) to obtain the displacement field u in terms of the phases θ n , we used the result of Eq. (6). This follows from the properties of moment tensors constructed from lattice vector sets Q = {q (n) } N n=1 . The p-th order moment tensor constructed from Q is given by In two dimensions, for a parity-invariant lattice vector set that has a B-fold symmetry, Ref. [Chen and Orszag, 2011] showed that all p-th order moments vanish for odd p and are isotropic for p < B. Every isotropic rank 2 tensor is proportional to the identity tensor δ i j , so for a 2D lattice vector set having four-fold symmetry, such as the set of shortest reciprocal lattice vectors {q (n) } 4 n=1 of the square lattice, we have 4 n=1 q (n) i q (n) j ∝ δ i j (Figs. 3 and 5 in Ref. [Skogvoll et al., 2021a] show the reciprocal lattice vector sets discussed in this appendix). Taking the trace and using that the vectors have the same length |q (n) | = q 0 , we get 4 n=1 q (n) i q (n) j = 4q 2 0 δ i j . In general, for any 2D parity invariant lattice vector set {q (n) } N n=1 with a B-fold symmetry where B > 2, we have 2D: As mentioned, this holds for the 2D square lattice, but it also holds for the 2D hexagonal lattice. In fact, the six-fold symmetry of the hexagonal lattice ensures that also every fourth-order moment tensor is isotropic, which results in elastic properties of the 2D hexagonal PFC model being isotropic [Skogvoll et al., 2021a].
To show this identity for a 3D parity invariant vector set with cubic symmetry, we generalize the proof in Ref. [Chen and Orszag, 2011] to a particular case of a 3D vector set that is symmetric with respect to 90 • rotations around each coordinate axis, such as the set of shortest reciprocal lattice vectors {q (n) } N n=1 of bcc, fcc or simple cubic structures. Let v be an eigenvector of Q i j with eigenvalue λ, i.e., is also an eigenvector of Q i j with the same eigenvalue λ. Repeating for a rotation around the y-axis demonstrates that Q i j has only one eigenvalue λ, so that it must be proportional to the rank 2 identity tensor Q i j ∝ δ i j . Taking the trace and using that the vectors have the same length |q (n) | = q 0 , we find 3D: Appendix C. Time derivatives of the dislocation density tensor Appendix C.1. Delta-function form Consider a moving dislocation line C = {r (λ, t)} of points r(λ, t) parametrized by the time t and a dimensionless λ which can be taken to go from 0 to 1 without loss of generality. Keeping the labelling fixed through its time evolution, we get α i j (r, t) = b j 1 λ=0 δ (3) (r − r (λ, t))(∂ λ r i (λ, t))dλ. (C.1) Suppressing the dependence of r on λ and t, we get taking the time derivative of Eq. (2), .
Inserting the expression for the delta-function in terms of the dislocation density tensor δ (2) (η n ) = α ik D (n) i q (n) k /(2π|D (n) | 2 ) into Eq. (14), we get Equating J (α) l j and J l j at a point r on the dislocation line, where α i j = t i b j δ (2) (∆r ⊥ ) using b · q (n) = 2πs n , lmn t m b j v n δ (2) (∆r ⊥ ) = Appendix E. Amplitude decoupling The (complex) polynomial f s (see Eq. (23)) results from the amplitude expansion of the ψ 3 and ψ 4 terms in Eq. (17). It may be computed by substituting Eq. (19) into Eq. (17) and integrating over the unit cell, under the assumption of constant amplitudes [Goldenfeld et al., 2005, Athreya et al., 2006, Salvalaglio and Elder, 2022. It features terms reading L =1 η n , with L = 3, 4 and n for which the condition L =1 q (n ) = 0 is satisfied. By multiplying this condition by b and using Eq. (8) it then follows that L =1 s n = 0. (E.1) In the equation for the dislocation velocity, Eq. (16), the only contributing amplitudes are those for which s n 0. The condition (E.1) implies that at least one of the other amplitudes, {η m } m n , appearing in terms of f s containing η n , also has s m 0 and then vanishes at the corresponding defect. Thus, for a given amplitude η n with s n 0, the terms in ∂ f s ∂η * n always contain at least one vanishing amplitude. Eq. (25) then reduces to Eq. (26) at the defect as η n = 0 and ∂ f s ∂η * n = 0 there. Importantly, a full decoupling of the evolution equation for amplitudes which vanish at the defect is obtained.