Properties of Dislocation Drag from Phonon Wind at Ambient Conditions

It is well known that, under plastic deformation, dislocations are not only created but also move through the crystal, and their mobility is impeded by their interaction with the crystal structure. At high stress and temperature, this “drag” is dominated by phonon wind, i.e., phonons scattering off dislocations. Employing the semi-isotropic approach discussed in detail in a previous paper (J. Phys. Chem. Solids 2019, 124, 24–35), we discuss here the approximate functional dependence of dislocation drag B on dislocation velocity in various regimes between a few percent of transverse sound speed cT and cT (where cT is the effective average transverse sound speed of the polycrystal). In doing so, we find an effective functional form for dislocation drag B(v) for different slip systems and dislocation characters at fixed (room) temperature and low pressure.


Introduction
Many modern material strength models, for example those in Refs. [1][2][3][4][5][6][7], are based on dislocation dynamics. However, dislocation mobility, especially in the high temperature and high stress regime, is poorly understood theoretically. Moving dislocations experience a drag due to their interaction with the crystal structure, and this drag coefficient B determines the dislocation glide time between obstacles. The lack of a well-established functional form for B(v, T, . . .) has led many researchers to assume B to be a constant (or a constant over a simple "relativistic" factor) as a fist order approximation within their strength models. Thus, better insight into the true functional form of B could improve those models.
Different mechanisms dominate dislocation drag in different regimes. However, at temperatures comparable to or higher than the Debye temperature and at high stress (leading to dislocation velocities in the range 0.01 v/c S < 1, where c S denotes the lowest shear wave speed corresponding to the direction of dislocation glide), phonons scattering off dislocations (commonly referred to as "phonon wind") constitutes the dominating effect. The lower end of this range is known as the "viscous" regime where B(v) at given temperature and pressure is known to be roughly constant. However, with increasing stress and thus increasing dislocation velocity, B exhibits a non-linear velocity dependence. This is seen in numerous molecular dynamics (MD) simulations (see e.g., [8][9][10][11][12] and references therein), but also within the recent theoretical framework in Refs. [13,14].
In Ref. [13], the theory developed by Alshits and collaborators [15] is taken to the next level by including not only the full velocity dependence of B, but also longitudinal phonons (in addition to the dominating contribution of transverse phonons) as well as an anisotropic dislocation field and single crystal elastic constants. Hence, this model developed having polycrystals in mind, keeps the phonon spectrum isotropic (for simplicity), but dislocations are modeled according to the single crystal symmetry (bcc, fcc, hcp, etc.) in order to take into account their anisotropy to some extent. This "semi-isotropic" approach constitutes an intermediate step in an ongoing long-term endeavor to include all anisotropic effects and the true phonon spectrum (which is beyond the scope of the current work). Nonetheless, valuable insights were already gained, e.g., the non-trivial dependence of the drag coefficient on the dislocation character angle ϑ (between line sense and Burgers vector).
For now, the model is also restricted to the subsonic regime where v < c T . The question whether dislocations in metals can reach supersonic speeds is still under debate, although numerous MD simulations suggest it is possible [8,11,[16][17][18][19]; see also the recent discussion on interpreting those results in the context of line tension and dislocation shape [20]. For recent literature on supersonic dislocations, (see, e.g., [21][22][23] and references therein).
Here, our main goal is to highlight the effective functional dependence of B on the dislocation velocity within the theory of [13] (and its numerical implementation of Ref. [24]), and to explain how to derive simple analytic representations of B(v) that are amenable to subsequent use in applications (such as material strength models). We also present new results for metals and slip systems not presented in [13] (i.e., prismatic and pyramidal slip for hcp metals). Thus, the current work in a sense complements Ref. [13] and the theory developed there.

Phonon Wind in the Semi-Isotropic Approach
The drag coefficient B of a dislocation is defined as the proportionality coefficient of the force F = Bv needed to maintain dislocation velocity v. It is related to the dissipation D per unit length due to phonon scattering via D = Bv 2 , and takes the form [25,26] where q and q are the wave vectors of incoming and outgoing phonons, s and s label their polarizations (two transverse and one longitudinal), and q is the wave vector associated with the dislocation field in Fourier space, d kk . In contrast to the phonons (which are quantized and have discrete wave vectors determined from the perfect lattice), the dislocation is modeled as a classical field in the continuum limit. Assuming an infinitely long, straight dislocation, its only spatial dependence is within the plane perpendicular to the dislocation line. The sums over discrete phonon momenta can subsequently be approximated as integrals over the first Brillouin zone. ω q , ω q are the phonon frequencies presently depending linearly on the wave vector length in accordance with the isotropic Debye approximation, i.e., ω q = c s | q | where c s is the sound speed of a phonon with polarization s (either transverse c T or longitudinal c L ). Ω q = | q · v| is the energy transfer whenever a phonon scatters on the dislocation. n q denotes the equilibrium phonon distribution function n q = (exp(hω q /k B T) − 1) −1 , which controls the number of scattering events per unit time.h is Planck's constant, ρ is the material density, and the two Dirac delta functions in the second line of Equation (1) encode momentum and energy conservation within each scattering event. Γ finally represents the associated matrix element, or scattering probability. As such, it depends on the (anisotropic) dislocation displacement gradient field d kk , the (quantized, isotropic) phonons whose orthonormal polarization vectors are presently denoted by w qi := w i ( q, s), and a linear combination of second (SOEC) and third order elastic constants (TOEC) of the anisotropic single crystal grains of a polycrystal,Ã i j k ijk . For technical details on the theory, we refer to [13] as well as [14,15].

Steady State Dislocations and Slip Geometries
The displacement gradient field in the continuum limit and within the realm of linear elasticity of a dislocation moving at constant velocity, can be determined from solving the equations of motion (e.o.m.) and the (leading order) stress-strain relations known as Hooke's law: where we have introduced the notation u k,l := ∂ l u k for the gradient of the displacement field u k , andü j := ∂ 2 u j ∂t 2 for the time derivatives. For constant velocity v i , this system of equations can be rewritten asĈ ijkl u k,il = 0 with "effective" elastic constantsĈ ijkl := C ijkl − ρv i v l δ jk (see [27]).
Upon introducing perpendicular unit vectors m 0 and n 0 , which are normal to the sense vector t of the dislocation, i.e., t = m 0 × n 0 , the solution takes the form u j,k (r, φ) =ũ j,k (φ)/r whereũ j,k (φ) is a function of Burgers vector, m, n, andĈ ijkl [28] (p. 476): with the shorthand notation (ab) jk := a iĈijkl b l . Variables r and φ are polar coordinates in the plane spanned by m = m 0 (ϑ) cos φ + n 0 sin φ and n = n 0 cos φ − m 0 (ϑ) sin φ, where n 0 is the slip plane normal and m 0 (ϑ) is perpendicular to n 0 and t(ϑ) = 1 b b cos ϑ + b × n 0 sin ϑ . As such, m 0 (ϑ) depends on the dislocation character angle ϑ and is parallel to v. The important feature to note is that u j,k (φ) includes terms proportional to (nn) −1 and hence exhibits divergences whenever det(nn) = 0. This happens at certain combinations of polar angle φ and critical velocity | v c |. As shown in Ref. [20], critical velocities are typically close to (and sometimes equal to) the lowest shear wave speed associated with the direction of v in the single crystal. All dislocation displacement gradients computed with the present method are hence restricted to (constant) velocities v that are smaller than v c .
As noted in the previous section, dislocation field u j,k (r, φ) (more precisely its Fourier transform d jk ( q)) enters Γ within Equation (1), and thus the drag coefficient B depends quadratically on u j,k . For simplicity, we presently only consider perfect dislocations; incorporating more realistic models of the dislocation core as well as the effect of partial dislocations into the dislocation drag coefficient are beyond the scope of the present paper and we leave those considerations to future work. For recent advances on the theoretical modeling of dislocation cores (albeit disconnected from phonon wind theory), see [29][30][31][32][33] and references therein.
The slip systems we have considered here are: where a and c are the lattice constants given in Tables 1 and 2, and b and n 0 denote the Burgers vector and slip plane normals, respectively (see Refs. [13,20] for details). For the case of close-packed hexagonal (hcp) crystals, we assume the basal plane is normal to the third axis in Cartesian crystal coordinates. The three hcp slip systems we consider, basal, prismatic, and pyramidal slip, share the same Burgers vector but have different slip plane normals. All except for the bcc slip system above lead to expressions that are symmetric with respect to ϑ → −ϑ, and all slip systems are π-periodic.
In Tables 1 and 2, we list all input data that were used in the computation of the drag coefficient below. For the effective Lamé constants of the polycrystal, we have chosen to use the separate experimental values (where available) listed in those tables rather than analytically averaging over the single crystal values. The only exceptions are Mo and Zr due to lack of experimental data, and because the Voigt and Reuss bounds are very close to each other in those cases. In fact, for SOEC of cubic crystals, analytic averaging would be a viable avenue as well (assuming negligible texturing), but not so much for hcp and other crystals, see [34] and references therein. (The single crystal averages for the Lamé constants of cubic crystals agree well -within a few percent -with the experimental results listed in Table 1, with the exception of Ni whose averaged shear modulus is ∼ 11% higher than the measured value, and also Au whose averaged λ is ∼ 12% lower than the measured value.) Table 1. List of input data for cubic crystals used in the calculation of the drag coefficient; all elastic constants are given in units of GPa. The references we used to compile these data are: Ref. [35] (Section 12) (lattice parameters a and densities ρ), Refs. [36] (p. 10) and [37] (effective Lamé constants of the polycrystal except for Mo), Ref. [35] (Section 12) (single crystal SOEC and Zener anisotropy ratio A := 2c 44 /(c 11 − c 12 )), and Refs. [38][39][40][41][42][43] (TOEC). The Lamé constants of Mo (marked with *) are analytical averages of the single crystal SOEC (see, e.g., [34]). The conventions for the single crystal elastic constants are those of Brugger [44].

The Low Velocity Limit
In the limit of small velocity v, small meaning v c T and v c S (where c T is the effective polycrystalline transverse sound speed and c S is the lowest shear wave speed of the single crystal in wherev denotes the unit vector in the direction of v. Explicit numerical calculations for a number of metals show that the first order velocity correction has a negative coefficient C < 0. To understand why this is the case, we note that the dislocation field itself depends only on the square of its velocity and thus its Taylor expansion around small v has no linear term. Furthermore, since ω q = c s q and Γ scales as 1/ω q ω q , the drag coefficient depends on the sound speeds as 1/c 3 s c 2 s . Since c T ∼ c L /2, the largest contribution to B at low velocity v is due to the purely transverse branch (where both incoming and outgoing phonons are transverse), as already observed in earlier work [13,14,25]. In this case, it is convenient to introduce a dimensionless integration variable proportional to the ratio t ∝ | q |/| q|. The energy conserving delta function then restricts the integration range of this new variable t such that it shrinks with growing dislocation velocity v (see [13,14]). This is the dominating effect and the reason for negative C.

High Velocity Limit
Our use of an isotropic Debye phonon spectrum introduces the limitation v < c T on our present theory. Nonetheless, B does not diverge at v = c T : All divergences within B are inherited from the poles present in the dislocation field, as pointed out in Section 2.1 above. Indeed, those appear at critical velocities v c , which depend on the slip geometry, material constants, as well as the dislocation character ϑ. To determine the highest degree of divergence, we first recall the study done in Ref. [14] in the purely isotropic limit and only for the transverse phonon modes: There it is found that the highest degree of divergence of a dislocation field for pure edge is 1/(1 − β 2 T ) m with β T := v/c T and m = 1 at polar angle φ = 0 (or π), whereas the one for pure screw exhibited the milder divergence of m = 1/2. Within B, where the dislocation field enters quadratically and angles φ are integrated over, this leads to initial estimates for the degree of divergence of B of m = 3/2 for pure edge and m = 1/2 for pure screw. However, within the purely transverse branch, the kinematic terms in Γ additionally suppresses the degree of divergence by 1, ultimately leading to B ∼ 1/(1 − β 2 T ) m as β T → 1 with m = 1/2 for edge and finite B for screw dislocations.
In the more general semi-isotropic case considered here, this latter cancellation cannot occur because now we have divergences at v c (ϑ), whereas the kinematic terms in Γ coming from the phonons only know about c T and c L . Similarly, the cancellation leading to the milder divergence of the pure screw dislocation in the isotropic limit is indeed special to the strictly isotropic case: For an isotropic (3) yield the milder divergence noted above. Finally, one must also not forget that edge and screw dislocations decouple only in the isotropic limit, but not in general, which is why mixed dislocations cannot be represented as superpositions of edge and screw in "real" crystals.
To sum up: we presently expect the highest degree of divergence of the drag coefficient c ) m with m = 3/2 for arbitrary dislocation characters ϑ. Indeed, this expectation is confirmed by numerical results, where the asymptotic region cannot be well represented by fitting functions with m < 3/2.

Results and Their Effective Functional Form
Based on the analysis of the previous section, the simplest form of a fitting function for the drag coefficient at fixed temperature, pressure, and dislocation character angle which captures its velocity dependence in the small v as well as in the asymptotic regime v → v c is given by As illustrated in Figure 1 for the example of Ni at room temperature and ambient pressure for a number of dislocation character angles ranging from pure screw (ϑ = 0) to pure edge (ϑ = π/2), Equation (6) is perfectly sufficient in some cases. Corresponding fitting parameters (in units of µPa s) and critical velocities-all dependent on ϑ-are listed in the figure legends and titles. However, if B shows a stronger v dependence in the intermediate region, which is the case for a number of metals and slip systems, additional terms are required to improve the fits. Candidates for such additional terms include of course any polynomial x k with k ≥ 2 or subleading divergences (which are always present), e.g., (1 − x 2 ) −m with 0 < m < 3/2 or ln(1 − x 2 ). Our goal is to keep B simple and the number of fitting parameters small. Empirically, we found that adding only one additional term, (1 − x 2 ) −1/2 , greatly improves the fits in most cases where Equation (6) is insufficient.
Hence, better fits to the drag coefficient from phonon wind for a dislocation of fixed character angle ϑ are achieved using the function Once again, it depends on the velocity in ratio to the critical velocity v c (ϑ). Note that, since (nn) is a 3 × 3 matrix, one always has three solutions for det(nn) = 0, and each can be represented as v c (φ). The branch with the smallest value for v c will lead to a divergence in (nn) −1 first. However, that solution need not always lead to a divergent drag coefficient since kinematics restrict the range of polar angle φ. This happens for example for pure screw dislocations in fcc metals where B diverges at a larger critical velocity v c than the dislocation field itself.  = v c /c T , and β T = v/c T . The solid lines show the results of numerically evaluating B according to Equation (1) using the software in Ref. [24], see Ref. [13] for details on the method.
According values for v c corresponding to divergences in B, as well as the five fitting parameters C i for pure screw and edge dislocations, and for B averaged over all dislocation character angles, and each for the various metals computed (at room temperature and ambient pressure), are listed in Tables 3-6. (Averages were computed as mean values from B(ϑ) for 91 character angles, 0 ≤ ϑ ≤ π/2 with B(−ϑ) = B(ϑ), for fcc and hcp metals, and from 181 character angles, −π/2 < ϑ ≤ π/2, for bcc metals using [24].) Comparisons of these fits to the numerically computed results for B are shown in Figures 2-5 as a function of velocity over effective transverse sound speed of the polycrystal, β T = v/c T . Fits using Equation (7) can of course be derived for any other character angle ϑ. All numerical results presented here can be reproduced with the software in Ref. [24] developed by the present author.
With the exception of Ag, Au, and Cd, most metals shown here have B well below 0.04 mPas in the low velocity regime, for the most part due to lower values of their transverse sound speeds c T (see Tables 3 and 4). Table 3. List of critical velocities v c (m/s), and fitting function coefficients C i (µPa s) for some cubic crystals. The critical velocities are given in units of m/s as well as in ratio to c T . The fits are only valid up to 0.99c T and do not capture the asymptotic behavior for those metals/dislocations whose critical velocity is larger than this value. Superscripts "e", "s", and "av" refer to "edge", "screw", and "average", respectively. Furthermore, v av c coincides with the smallest critical velocity for all dislocation characters ϑ within the slip system.  Table 4. List of critical velocities v c (m/s), and fitting function coefficients C i (µPa s) for basal slip. The critical velocities are given in units of m/s as well as in ratio to c T . The fits are only valid up to 0.99c T and do not capture the asymptotic behavior for those metals/dislocations whose critical velocity is larger than this value. Superscripts "e", "s", and "av" refer to "edge", "screw", and "average", respectively. Furthermore, v av c coincides with the smallest critical velocity for all dislocation characters ϑ within the basal slip system.  Table 5. List of critical velocities v c (m/s), and fitting function coefficients C i (µPa s) for prismatic slip. The critical velocities are given in units of m/s as well as in ratio to c T . The fits are only valid up to 0.99c T and do not capture the asymptotic behavior for those metals/dislocations whose critical velocity is larger than this value. Superscripts "e", "s", and "av" refer to "edge", "screw", and "average", respectively. Furthermore, v av c coincides with the smallest critical velocity for all dislocation characters ϑ within the prismatic slip system.  Table 6. List of critical velocities v c (m/s), and fitting function coefficients C i (µPa s) for pyramidal slip. The critical velocities are given in units of m/s as well as in ratio to c T . The fits are only valid up to 0.99c T and do not capture the asymptotic behavior for those metals/dislocations whose critical velocity is larger than this value. Superscripts "e", "s", and "av" refer to "edge", "screw", and "average", respectively. Furthermore, v av c coincides with the smallest critical velocity for all dislocation characters ϑ within the pyramidal slip system.  (6)), especially in light of the uncertainty in our current model for B, which is hard to quantify. (In our example of nickel, both Equations (6) and (7), yield exactly the same fit for pure edge dislocations, whereas, in the case of screw dislocations, the four-parameter Equation (7) slightly improves an already decent fit by making use of the additional term 1/ √ 1 − x 2 , thereby changing also the values of the three other fitting parameters.)

Pyramidal
For one, we considered only isotropic phonons, whose spectrum deviates from the true one especially in the high frequency regime. Furthermore, we have neglected the dislocation core as well as the separation of dislocations into partials. However, the uncertainties in the experimental (or computational) determination of the TOEC also have a large effect on the accuracy of our present predictions. We also need to stress that the present model is limited to the subsonic regime, i.e., v < v c (ϑ) and v < c T . Furthermore, only straight dislocations moving at constant velocity were considered, i.e., the effect of acceleration or changes in shape are not (yet) considered. Finally, the stress field required to reach velocities close to v c will likely lead to sizeable temperature and pressure gradients, which would have to be considered in future improvements to B as well. A first attempt at incorporating the temperature dependence into B is currently in progress [50]. Direct comparison of B to experiments is limited to the low velocity regime (low meaning the viscous regime of β T ∼ 0.01): As (in part) pointed out in Ref. [13], our predictions for B(β T ∼ 0.01) agree well with experimental results for Al (ranging from ∼ 0.005 mPas to ∼ 0.06 mPas, cf. [51][52][53]) and Cu (ranging from ∼0.0079 mPas to ∼0.08 mPas, cf. [54][55][56][57][58]. MD simulation results are in the range ∼0.007-0.2 mPas for Al [8,59,60], and ∼0.016-0.022 mPas for Cu [10,61]. Our predictions are lower than experimental results for Fe (∼0.34 mPas for edge and ∼0.661 mPas for screw, cf. [62]) and Zn for both basal slip (0.035 mPas for edge and ∼0.034 mPas for screw, cf. [63]) as well as for pyramidal slip (0.27 mPas for edge and ∼0.16 mPas for screw, cf. [64]).
Our drag coefficient for Mo is lower than the MD-simulation value of ∼0.078 mPas for edge dislocations reported in [65]. Similarly, our drag coefficient for Ni is lower than the MD-simulation results of 0.0321 mPas for edge dislocations reported in [65], and ∼0.015 mPas for edge dislocations reported in [8,12,66], albeit close to the latter.

Conflicts of Interest:
The author declared no conflict of interest.