Properties of Accelerating Edge Dislocations in Arbitrary Slip Systems with Reflection Symmetry

We discuss the theoretical solution to the differential equations governing accelerating edge dislocations in anisotropic crystals. This is an important prerequisite to understanding high-speed dislocation motion, including an open question about the existence of transonic dislocation speeds, and subsequently high-rate plastic deformation in metals and other crystals.


Introduction and Background
Dislocations can influence many materials' properties such as thermal conductivity [1], thermal stability [2], impact sensitivity [3], ferroelectricity [4], and electrical resistance [5]. At extremely high rates, plastic deformation is governed by high-speed dislocations, a regime where dislocation mobility is poorly understood [6][7][8]. High-speed dislocations experience a drag force due to scattering phonons (an effect known as 'phonon wind') and this interaction (and thus, dislocation mobility) is sensitive to the stress distribution in the vicinity of the moving dislocation. Dislocation drag is thus key to predicting material strength at extremely high stress and understanding high-rate plastic deformation [9]. The first-principles phonon wind theory was derived in the isotropic and steady-state limit for dislocation glide velocities that are much smaller than the transverse sound speed some time ago; see the excellent review article [10]. More recently, dislocation drag theory was generalized to very high (but still subsonic) dislocation velocities [11] and anisotropic crystals [12], though the effects of acceleration have so far been neglected.
Another key question in this regard is whether dislocations can reach transonic and supersonic speeds under sufficiently high stress. The only indication that such speeds are possible comes from molecular dynamics (MD) simulations [13][14][15][16][17][18][19]. Experiments cannot track dislocations in real time at these high speeds (After the present manuscript was completed, ref. [20] appeared, which for the first time measured transonic dislocations in diamond in real time), but one can hope to indirectly determine the presence of supersonic dislocations and perhaps estimate the fraction and velocity of these dislocations in the near future [21,22]. This in turn requires a thorough understanding of the solutions to the differential equations governing dislocations, i.e., the equations of motion supplemented by the (leading order) stress-strain relations. Dislocation theory predicts divergences in self energy and stress at certain limiting velocities [23][24][25][26] for steady-state dislocations. In the isotropic limit, it was shown [27][28][29] that an acceleration term together with a regularized dislocation core removes the divergence, thereby opening the possibility of supersonic events. Other authors emphasized the importance of size variations of the dislocation core as a function of dislocation velocity [30][31][32][33]. The steady-state solution for dislocations in arbitrary anisotropic crystals has been known for some time [34,35]. The case of accelerating dislocations in anisotropic crystals has also been studied [36][37][38][39][40][41][42], with pure screw dislocations having been discussed in the most detail [36,42]. The most general solution has been given only in a very formal form [39], apart from an additional asymptotic wave front analysis. In this paper, we consider formal derivation of the accelerating dislocation field of ref. [39] as a starting point to discuss in detail the solution of an accelerating pure edge dislocation in anisotropic crystals.
In particular, we discuss the solution to the following set of differential equations for accelerating dislocations for the special case of pure edge dislocations: in coordinates aligned with the dislocations, i.e.,ẑ is aligned with the dislocation line andŷ is parallel to the slip plane normal. The components of the tensor of second-order elastic constants (SOEC) is always measured in Cartesian coordinates that are aligned with the crystal axes, and thus this tensor must be rotated into our present coordinate basis, i.e.: with rotation matrix U.
In order to study pure edge (or pure screw) dislocations, the rotated tensor of SOEC must fulfill the following symmetry requirements (shown here in Voigt notation which maps index pairs to single digits, (11, 22, 33, 32/23, 31/13, 21/12) → (1, 2, 3, 4, 5, 6)): i.e., the six components c 14 , c 15 , c 24 , c 25 , c 46 , and c 56 must vanish, see refs. [43,44]. This ensures that u 3 = 0 implies ∂ i σ i3 = 0, and likewise that u so that there exists a u 3 that solves the differential equations independently from the pair (u 1 , u 2 ) and vice versa. Note that in the present coordinates, u i can only depend on x, y, and t, but not on z. This latter property implies that non-vanishing components c 34 and c 35 are allowed since they do not enter the differential equations above for pure screw or pure edge dislocations. On the other hand, the stronger condition c 34 = 0 = c 35 implies that the x 1 , x 2 plane is a reflection plane (and then σ 33 = 0 for pure screw dislocations rather than the weaker ∂ 3 σ 33 = 0). The most general solution for pure screw dislocations was recently derived in ref. [42]. The case of accelerating pure edge dislocations was previously studied by Markenscoff and Ni for the special case of c 16 = 0 = c 26 (in addition to (3)) in refs. [37,38], and the general case was presented in ref. [39]. In refs. [37,39], only a formal solution was derived, though not in closed form. Here, we present for the first time, a numerical implementation of the accelerating dislocation field for pure edge dislocations in various anisotropic slip systems and study its properties. Our code is included in version 1.2.7 of PyDislocDyn [45].

Most General Differential Equations for Pure Edge Dislocations
Following ref. [39] in this subsection, but setting u 3 = 0 and plugging the most general rotated tensor of SOEC fulfilling the required properties for studying pure edge dislocations, Equation (3), into the differential Equation (1), we find: Note that we have dropped the primes on the elastic constants for notational simplicity; nonetheless, all c ij are understood to be in the rotated frame aligned with the edge dislocation. Additionally, we have the boundary conditions where Θ(x) denotes the Heaviside step function, b is the Burgers vector length, and the slip plane is located at y = 0. Clearly, the above differential equations and their boundary conditions simplify significantly when c 16 = 0 = c 26 , which is what was studied in ref. [37,38].
In order to solve these more general equations, we apply a Laplace transform in time, i.e., as well as a two-sided Laplace transform (which is related to the Fourier transform with sλ → ik) in x, i.e., and thus, U i (λ, y, s) ≡ T {L{u i (x, y, t)}}. The transformed differential equations read Likewise, the transformed boundary conditions in the upper half plane (y ≥ 0) read where η(x) ≡ l −1 (x) and the integral over time was carried out explicitly as described in ref. [42]. Additionally, we demand lim [39] argues that the problem can be reduced to a problem on a half-space, so that we now assume y ≥ 0 in the following derivation, and we will generalize to negative y only at the very end. Note that the first term in boundary condition (9) is identified as that of the static problem which cannot be treated by a Laplace transform without running into convergence issues [42,46]. Hence, we presently subtract the static contribution and will add it at the end of our derivation; more precisely, we will add the well-known solution to the static problem at the very end so as not to clutter our equations in intermediate steps. Focusing only on the dynamic part of the accelerating dislocation field, we presently replace (9) with and for notational simplicity, we drop the tilde below ( U 1 → U 1 ). We furthermore assume that c 12 + c 66 = 0, i.e., we do not include the so-called irregular hyperbolic case [38] in our discussion, as we are unaware of any slip systems that in practice would exhibit this property [23]. The differential Equation (8) can be rewritten in 4 × 4 matrix form as where we defined the compliances as s ki C i22j ≡ δ kj , i.e., Since we focus here on the regular hyperbolic case, we may assume that the eigenvalues of the so-defined 4 × 4 matrix (µ m with m = ±1, ±2) are distinct [39]. Given these eigenvalues, we make the ansatz Plugging this ansatz into the differential Equation (8) yields the determinantal equation which may be used to calculate the µ m (λ) by solving the following fourth-order polynomial: Note that s is factored out in this equation so that µ m depends on λ but not on s. Finally, the asymptotic condition lim y→∞ ∂ 2 U i = 0 tells us that the sum over m in the ansatz (14) above must only include the positive eigenvalues and Markenscoff argued in [39] that because the slowness surface (whose equation coincides with the determinantal Equation (15) above) is symmetric about the origin, there are presently two positive eigenvalues, m = 1, 2. The corresponding eigenvectors are (A 1m , A 2m , −µ m sA 1m , −µ m sA 2m ) where the A im is determined from together with the boundary conditions which presently read Plugging the ansatz A 2m = a m A 1m into (17), we find for a m : where the last equality follows from the fact that µ m solves (15). The boundary conditions (18) finally determine A 1m , and written in matrix form we presently have Thus, with a m given in (19). Note that the coefficients A im (λ) do not depend on s; this will be important later when we derive the inverse Laplace transform.

Cagniard-De Hoop Method
In order to determine the displacement gradient field in real space and time, we need to apply the inverse Laplace transform +i∞ −i∞ f (λ)e −sλx sdλ and integrate λ along a line parallel to the imaginary axis. This latter integral will not be carried out explicitly, but rather we want to rewrite it in a way that allows us to interpret this integral as a Laplace transform in time so that a subsequent inversion of the one-sided Laplace transform L{u i } need not be carried out explicitly.
Thus, for each term in U i we interpret the following combination as a strictly positive time variable τ in order to apply the Cagniard-De Hoop method [47][48][49]: The reader is reminded that we presently restrict our calculation to the half plane y ≥ 0. In order to be able to integrate τ over the positive real axis instead of over the imaginary λ axis, one needs to study an integral over λ over a closed path in complex space and to account for the residua of all enclosed poles. This step requires knowledge of the locations of all poles in the expressions above, and hence knowledge of the roots µ m (λ). Note, that such poles occur only for transonic and supersonic dislocations, but not in the subsonic regime [39]. Furthermore, in passing from integration variable λ to integration variable τ m , we need the inverse of function (22), i.e., λ m (τ m ), as well as the Jacobian dλ m dτ m . The inverted functions λ m appear in complex conjugate pairs which both need to be taken into account in order to integrate over a closed path [39,42]. Using Cauchy's theorem, we presently have in the subsonic regime: where τ min m = lim λ→0 τ m (λ) and A jm =Ã jm U 0 is given in (21) with (19). In the transonic and supersonic regimes, the expression above needs to be supplemented by appropriate residua for all enclosed poles in the integration path. As discussed in earlier papers [42,46], calculating u j directly is troublesome due to subtleties with respect to poles, and it is generally better to solve for its gradient. Thus, taking derivatives with respect to x and y prior to passing from λ to τ, we find Another important subtlety concerns the exchange of integrals over λ and x prior to the change of variables, which is only permissible if both integrations converge absolutely; this is not the case in general and a remedy was put forward in the context of pure screw dislocations in refs. [42,46]. In particular, the exchange of integrals leads to poles on the slip plane at y → 0 which stem from the first two terms of a Taylor expansion of η(x ) around x = x. On the other hand, if one were to replace η with its linear order Taylor expansion terms, the integral over x can be carried out analytically before changing integration variables: In that case, τ will not depend on x (i.e., one defines (22) with x = 0) and only one integral over λ (resp. τ m ) is left.
To sum up: In order to eliminate divergences on the slip plane in the x integration, we must add and subtract the dynamic term with η(x ) replaced by its linear order Taylor [42]. Hence, Considering the properties of the Laplace transform, where multiplication by e −sT corresponds to a translation in time t → t − T and multiplication by s corresponds to a time derivative (modulo boundary terms which are zero here), we can read off the solution: where λ m depends on the appropriately shifted time , matching in each term the according part of the argument of the step function.

Special Cases: Constant Velocity and Constant Acceleration Rate
The simplest case one can study within the present solution is a dislocation initially at rest which suddenly starts moving at constant velocity v at time t ≥ 0. As discussed previously in the context of pure screw dislocations in [42], this "jump" in velocity is nonphysical, but in the large time limit the solution must tend to the well-known steadystate solution, thus providing us with a consistency check. The assumption of constant dislocation velocity at t ≥ 0 leads to the following simplifications: Due to the last equality, the second and fourth lines within Equation (27) (i.e., the terms containing the time derivative and the integral over x ) vanish identically for a dislocation moving at constant velocity.
The simplest physical case within the present dynamic solution follows from the assumption that the dislocation is at rest at time t < 0 and starts to accelerate at a constant rate a from time t ≥ 0. Then l(t) = a 2 t 2 > 0 and hence [42] η The velocity at time t is given by v(t) = at and the transition from subsonic to transonic happens when t = v lim /a, where v lim is the lowest limiting velocity whose value can easily be computed using the review article [23] and/or the open source code [45].
We have implemented this constant acceleration rate case in Python, using a combination of symbolic (sympy) calculations and numerical methods, and have integrated it into the code PyDislocDyn [45]. The general strategy is as follows: The material's tensor of elastic constants is rotated into coordinates where the dislocation line is parallel to the z direction, the slip plane normal points in the y direction, and the edge dislocation accelerates from rest in the x direction at rate a. We then calculate the time t 1 at which the accelerating dislocation reaches a user-specified target velocity, as well as the position of the dislocation core at that time in order to shift the x coordinate such that the dislocation core resides at the origin at time t 1 . We use sympy to calculate the four solutions µ(ρ/λ 2 ) from Equation (16) after plugging in numerical values for all (rotated) elastic constants and the material density, i.e., λ is the only unknown. For each of these 4 solutions, we determine τ(λ) and its derivative, and the resulting sympy expressions are subsequently 'lambdified', i.e., converted into functions of λ, x, and y. We then loop over all points x, y we wish to determine the displacement gradient for. At a given point x, y, function τ depends only on λ, and since we are interested in one snapshot in time (meaning we know τ), we can numerically determine λ(τ); note that λ is a complex number and we use mpmath's recommended root-finding method (the Muller method). This step constitutes the bottleneck of our implementation, i.e., calculating the dislocation field for accelerating edge dislocations is orders of magnitude slower than for screw dislocations which were discussed in [42]. Once we have λ, we determine µ(λ) and the Jacobian 1/ dτ dλ . At this point, we have four sets of λ, µ(λ), but only two satisfy the asymptotic condition lim y→∞ ∂ 2 U i = 0. Markenscoff [39] determined that the imaginary parts of λ and µ/λ must have opposite signs for positive y, and we drop the other two solutions to λ. The remaining two sets of λ, µ(λ) are plugged into (19) and (21), and subsequently into the first (i.e., leading) dynamic terms of (27). The static part is computed with the well-known Stroh/integral method [34]. The time-derivative term in (27) can be neglected for constant acceleration rates. Figure 1 shows the edge dislocation field at the example of hcp Mg for prismatic slip and compares the accelerating field to the steady-state field. Figure 2 shows the edge dislocation field at the example of bcc Nb for the 112 slip planes and compares the accelerating field to the steady-state field. In contrast to the previous example, edge dislocations on 112 slip planes of bcc metals have a non-vanishing (rotated) elastic constant c 26 , and thus represent a more general case than the former. Both examples show some enhancement of the dislocation displacement gradient field for moderate acceleration rates of a ∼ 10 13 m/s 2 typical for flyer plate impact scenarios [50], albeit maintaining the shape of the steady-state solution for the most part. Only for very extreme acceleration rates do we start to see more significant deviations as illustrated in Figure 3 with the example of Mg. Note that the numerical accuracy of the accelerating edge solution in its current implementation is limited by the accuracy of the (complex) root-finding algorithm.  [51]). This velocity corresponds to roughly 92% of the critical velocity. All plots are centered at the dislocation core, showing the plane perpendicular to the dislocation line in units of a Burgers vector. On the left of each pair of plots, we show the steady-state solution [34] and on the right we show the full solution for constant acceleration (27) with (29) and a = 1 × 10 13 m/s 2 at time t v = v/a = 2.838 × 10 −10 s needed to reach velocity v. At this point, the dislocation has traveled a distance of 0.4 microns. We see that the changes in the dislocation displacement gradient due to the inclusion of acceleration lead to a slight enhancement.  [51]). This velocity corresponds to roughly 90% of the critical velocity. All plots are centered at the dislocation core, showing the plane perpendicular to the dislocation line in units of a Burgers vector. We compare the steady-state solution [34] with the full solution for constant acceleration (27) with (29) and a = 1 × 10 13 m/s 2 at time t v = v/a = 1.883 × 10 −10 s needed to reach velocity v. At this point, the dislocation has traveled a distance of ∼0.18 microns. We see that the changes in the dislocation displacement gradient due to the inclusion of acceleration lead to a slight enhancement. Furthermore, we confirm (numerically) that the divergence at a 'critical' dislocation velocity (which separates the subsonic from the transonic regime) persists for general accelerating edge dislocations with vanishing core size, consistent with previous work on the isotropic limit [27] as well as the accelerating screw dislocation in anisotropic crystals [42].

The Isotropic Limit
The following simplifications apply in the isotropic limit: c 22 = c 11 = c 12 + 2c 44 , c 66 = c 44 , and c 16 = 0 = c 26 , as well as s 11 = 1/c 11 = s 22 and s 12 = 0 = s 21 within (13). Hence, Equation (16) simplifies to where c 11 = c 12 + 2c 44 , and solutions µ m are found to be In both cases, only one of the two signs must be considered, namely, convergence of (14) requires that the real part of µ m has the same sign as y. For positive y, this means that Im(λ) > 0 implies Im(µ m /λ) < 0 and vice versa [39].
The definition of τ m (with x = 0) then yields with r 2 ≡ x 2 + y 2 and the short-hand notation c 1 ≡ c T = c 44 /ρ and c 2 ≡ c L = c 11 /ρ for the transverse (T) and longitudinal (L) sound speeds. This special case was discussed in ref. [52]. If we assume a constant dislocation velocity from time t > 0, i.e., η(x) = x/v and take the limit of t → ∞ after translating our coordinates to move with the dislocation (i.e., replacing x = x + vt, r 2 = (x + vt) 2 + y 2 everywhere prior to taking the limit, see [42]), we recover the well-known steady-state solution for an edge dislocation in an isotropic medium [11,53]:

Conclusions
In this paper, we have presented and discussed the full solution to the differential equations for an accelerating edge dislocation in a general anisotropic crystal in the subsonic regime. Taking the formal solution of ref. [39] one step further, we have derived the edge dislocation displacement gradient field using a combination of analytical and numerical methods. Our python implementation is included in version 1.2.7 of the code PyDislocDyn [45]. Two examples were illustrated in Figures 1 and 2 showing that the dislocation strain field is slightly enhanced in the accelerating case, at least for typical dislocation acceleration rates of a ∼ 10 13 m/s 2 [50], though still similar enough to the steady-state solution (except for extreme conditions such as very high acceleration rates and velocities near the limiting velocity), so that in most larger simulations it makes more sense to use the (several orders of magnitude) faster-to-compute steady-state solution. The transonic regime of the accelerating edge dislocation as well as accelerating mixed dislocations are left for future work.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.