Statistical mechanics of kinks on a gliding screw dislocation

The ability of a body-centered cubic metal to deform plastically is limited by the thermally activated glide motion of screw dislocations, which are line defects with a mobility exhibiting complex dependence on temperature, stress, and dislocation segment length. We derive an analytical expression for the velocity of dislocation glide, based on a statistical mechanics argument, and identify an apparent phase transition marked by a critical temperature above which the activation energy for glide effectively halves, changing from the formation energy of a double kink to that of a single kink. The analysis is in quantitative agreement with direct kinetic Monte Carlo simulations.

The rate of plastic deformation of a body-centred-cubic (bcc) metal depends strongly on temperature, with pure iron notoriously becoming brittle below freezing [1]. The temperature of this brittle-to-ductile transition, where the metal is brittle at a low temperature and is ductile at an elevated temperature, is raised by up to hundreds of degrees after irradiation by highly energetic particles, which leads to stringent requirements on the minimum operating temperature of ferritic steels and other bcc materials for technological applications in a radiation environment [2][3][4]. In order to predict how the mechanical properties vary over the service lifetime of a structural component exposed to radiation at elevated temperatures, it is desirable to develop an explicit model, relating microstructural changes to plastic deformation.
The principal rate-limiting mechanism for plastic deformation in bcc metals is the thermally activated motion (glide) of screw dislocations [5], which are topological line defects in the crystal lattice, acting as carriers of plastic deformation. A screw dislocation can be approximated by an elastic line placed in a periodic Peierls potential [6,7], driven by shear stress and by random thermal noise from the surrounding crystal lattice [8,9]. At low stresses, screw dislocations move by the kink mechanism, in which thermal fluctuations first lead to the nucleation of a pair of kinks into the next valley of the Peierls potential, which then subsequently diffuse athermally along the dislocation [10] to advance the rest of the line by a discrete amount. As a result, one would expect the activation energy characterising the plastic strain rate to be close to the formation free energy of a kink pair 2f k . However, in an apparent contradiction to this argument, * max.boleininger@ukaea.uk † gallauer@maths.ox.ac.uk ‡ sergei.dudarev@ukaea.uk § swinburne@cinam.univ-mrs.fr ¶ daniel.mason@ukaea.uk * * danny perez@lanl.gov . Right: Entropy as a function of reciprocal temperature β = (kBT ) −1 computed using the free energy model (19). The entropy becomes extrinsic for temperatures above a certain critical value, indicated by a dashed line in the right panel, see equation (8).
experimental observations of the brittle-to-ductile transition in high-purity bcc metals appear to exhibit activation energies closer to the formation free energy of a single kink f k [1,[11][12][13]. Swinburne and Dudarev [13] offered a resolution to this contradiction by proposing a critical transition in the activation energy for dislocation glide. By analyzing kink populations on a dislocation in thermal equilibrium, they concluded that a dislocation segment longer than a certain temperature-dependent critical length L * glides with an activation energy of f k , while shorter dislocation segments glide with an activation energy of 2f k . This conclusion is also reached in the kink diffusion model by Hirth and Lothe [14], based on an argument that kinks start annihilating if the line length exceeds the mean spacing between thermal kinks ∼ L * [15,16]. The segment length-dependent mobility affects the ability of a screw dislocation to unpin from obstacles [13], and as such has major implications for the predictive interpretation of experiments on obstacle hardening and embrittlement. Yet despite the extensive use of screw dislocation mobility laws in coarse-grained methods for modelling plastic deformation [17][18][19][20], to our knowledge there is still no suitable analytical expression able to capture the full complexity of dislocation mobility, consistent with the microscopic statistical mechanics description of a fluctuating dislocation line, including the critical transition noted above.
In this letter, we give a closed-form expression for the glide velocity of a screw dislocation in the kinklimited regime as a function of temperature, resolved shear stress, and segment length. This analytical model explicitly captures the critical transition, and is in quantitative and qualitative agreement with kinetic Monte Carlo (kMC) dislocation dynamics simulations.
The model describes a screw dislocation, moving in a glide plane, by a set of n ∈ Z + displacement variables u i adopting integer values and satisfying the Born-von Karman periodic boundary condition u n+1 = u 1 . The dislocation line is spanned by points (ib, hu i ), where i is an integer, b is the Burgers vector length and h is the distance between the adjacent Peierls valleys. A site i may accommodate any number of left kinks (u i+1 − u i > 0) or right kinks (u i+1 − u i < 0), or no kinks at all (u i+1 − u i = 0). The positioning of kinks along the line is enough to uniquely characterise a line configuration, and in what follows we define the set of microstates as the set of all the unique line configurations. Some representative line configurations are illustrated in Fig. 1.
The coarse-grained Hamiltonian of a dislocation line in the absence of external stress is where f k (T ) is half the free energy of formation of a kink pair [21] and r is the number of kink pairs on the line. In this representation a kink is described as a quasiparticle with a distinct free energy of formation f k (T ), with the dependence on temperature originating from the coarsegrained atomic degrees of freedom [22]. This model is closely related to the discrete Gaussian model [23]. Equilibrium properties.-The atomic lattice acts as a heat bath with temperature T . In the absence of applied stress, biasing the motion of the dislocation line, or in other words in the limit σ = 0, the microstates containing identical numbers of kinks are degenerate. The system equilibrates, with thermal fluctuations giving rise to a balanced rate of creation and annihilation of paired left and right kinks. Placing first r left kinks on n available sites, and subsequently r right kinks on n − r remaining sites, we find the canonical partition function where z = exp(−2βf k ), β = (k B T ) −1 is the reciprocal temperature, k B is the Boltzmann constant, and x is the floor function. The above summation result is exact, and is expressed in terms of a hypergeometric function 2 F 1 (a, b; c; x) [24]. A similar partition function was previously investigated in the limiting case where left and right kinks were assumed to be indistinguishable [13]. The exact solution expressed as a hypergeometric function is not conducive to further analysis. Hence in what follows we restrict the discussion to the kink-dominated mobility regime, where z 1. In this limit, pertinent to virtually all the conditions encountered in experiment [25], the partition function is given by the expression where I 0 (x) is the modified Bessel function of the first kind of order zero. We refer to Appendix A for the detailed derivation.
Since the system studied here has a finite size n, the entropy S = ∂ T (k B T ln Z) is not necessarily extrinsic. This is particularly evident in the low temperature (LT) limit, where In the high temperature (HT) limit, entropy has the form As a function of system size and temperature, entropy changes from being intrinsic (S LT ∝ n 2 z) at low temperature to being extrinsic (S HT ∝ n √ z) at sufficiently high temperature. Furthermore, the argument of the exponential function, dominating the variation of entropy as a function of temperature, effectively halves in the high T limit, from 2f k (where S ∼ z = exp(−2βf k )) at low temperature to f k (where S ∼ √ z = exp(−βf k )) at high temperature, in agreement with Refs. [13,14,26].
In what follows, we define the critical temperature T * as a temperature of a transition between the intrinsic and extrinsic regimes for a given dislocation segment size n. To illustrate the nature of this transition, consider for example the mean number of kink pairs r on the dislocation line where I 1 (x) is the modified Bessel function of the first kind of order one. Similarly to entropy, the expression for ln r exhibits a change of slope as a function of reciprocal temperature. β * = (k B T * ) −1 corresponds to the point where the change of slope is maximum: The value of T * defined by this equation depends both on the length of the dislocation segment n and the free energy of formation of a kink pair 2f k . The resulting expression for the root of equation (7) does not admit an analytical solution as it depends on the free energy function f k (T ). Approximating the free energy of formation of a kink pair by a linear function of temperature f k (T ) ≈ a 0 + a 1 T , which is motivated by atomistic free energy simulations in tungsten [22], we find that the critical temperature satisfies the following implicit equation where ξ = 0.95483 is a numerical constant. The Arrhenius plot in Fig. 2 shows the dislocation segment size scaling of the transition. We note that this critical transition is not a conventional thermodynamic transition, as it occurs in a finite size system. A related apparent roughening transition is also found in the sine-Gordon model [27], which is a continuum equivalent of the discrete Gaussian model [28].
Non-equilibrium steady state.-The introduction of bias in the form of applied shear stress drives the dislocation line out of equilibrium. The rate of nucleation and propagation of kinks in the direction favoured by the applied bias is higher, while the processes occurring against the bias are suppressed. The presence of bias eventually gives rise to the steady-state drift of the dislocation line. With no loss of generality we restrict the following discussion to the case σ > 0.
To evaluate the mean glide velocity, one needs to know the probability of occurrence of any microstate and its associated rate of escape into the connected microstates. We describe the kinetics of the line as a sequence of fundamental reactions, where each involves moving just one line segment at a time, either in the positive (u i → u i +1) or negative (u i → u i − 1) sense with respect to the biasing direction. Each positive or negative reaction thereby changes the mean position of the line by ±h/n. The reaction rates are distinguished by whether the reaction results in motion of a kink or creation of a kink pair: 1. If the reaction u i → u i ± 1 increases the number of kinks on the line, a kink pair is created with rate k ± c . 2. If the reaction u i → u i ± 1 conserves or reduces the number of kinks on the line, a kink is moved with a rate k ± m . For the sake of simplicity we assumed that the process of a kink pair annihilation is equivalent to a right kink moving to a site containing a left kink, or the other way around.
We begin by expressing the non-equilibrium dynamics of the system in terms of a Markov chain [29], as this will allow to systematically introduce the approximations required to arrive at an analytical expression for drift velocity. Let the state vector (π m ) i contain the relative number of visits to microstate C i after m reactions. The normalization i (π m ) i = 1 applies because no sinks or sources are present. Depending on how the kinks are distributed in C i , a microstate may reach a connected microstate C j through one of the fundamental reactions with rate k ij . Defining k i = j k ij as the total escape rate out of C i , the transition matrix is obtained as (P) ij = k ij /k j . The state vector is then iterated according to: The state vector converges independently of the initial value π 0 to the same steady state lim m→∞ π m = π. Expectation values x in the steady state are obtained by the mean outcome x i weighted by the average time π i k −1 i spent in a given microstate : where p i is the probability for the system to be in microstate C i . The outcome for velocity v i of a microstate C i is obtained as the net balance between escape rates associated with reactions moving the line in the positive or negative sense with respect to the direction of applied stress where n ± i,ε is the number of microstates connected to C i that can be reached by a fundamental reaction occurring at rate k ± ε . The transition matrix can only be built explicitly for small systems as the number of microstates scales exponentially with n. An exact solution for the drift velocity is consequently unfeasible. Instead we introduce a series of approximations that eventually lead to an analytical solution.
First, we reduce the dimension of the state vector by partitioning microstates with similar probability of occurrence into non-overlapping groups P α , each containing a certain number of microstates N α . Next, we assume that the partitioning was chosen such that all the microstates C i ∈ P α in a group occur with approximately equal probability p i ≈ N −1 α p α , namely where is the group-averaged outcome.
In practice, we choose to put together all the microstates with the same number of kink pairs r into a group P r . This is motivated by the fact that the rate of kink motion is much higher than the rate of kink creation, allowing the system to rapidly explore microstates with the same kink number between the reactions leading to the creation of kink pairs. The description could be improved by further grouping the states by area A, but this results in a severe combinatorial challenge [30]. In the limit where the bias is relatively weak, the occupation probability p r of a group P r will be very close to the equilibrium occupation probability p eq r . Equating the identity with z replaced by the normalized average rate of creation of kink pairsz = (k + c + k − c )/(2k 0 ) to approximately account for the biased creation of kink pairs in the preferred direction, with k 0 = k ± c | σ,T =0 . Note that the linear-response formalism is recovered by choosing z = z. In practice we find that the average rate gives an improved agreement at high stresses, see Appendix for a comparison.
Following Eq. (13), we proceed to determine the groupaveraged line velocity v r . To achieve this, we need to evaluate the total number of configurationally permitted escape reactions from a set of microstates in P r . The resulting combinatorial problem can be solved exactly, provided that we restrict a site on the line to accommodate at most one kink. The exact solution is possible, it involves hypergeometric functions and is presented in the Appendix for the sake of completeness. Here we continue with one final approximation. Bearing in mind that the approximations taken so far are valid in the kinkdominated regime z 1, where the kink density is low [25], the number of microstates where kinks overlap is very small in comparison with the number of microstates where kinks do not overlap. Hence the number of microstates connected to a given group P r can be well approximated by considering every site to be free for kink nucleation (n ± r,c ≈ n), and every kink to be free to move in either direction along the line (n ± r,m ≈ 2r). Using this argument, we arrive at an expression for the group-average velocity in the form where N r = n 2r /(r!) 2 is the number of microstates containing r kink pairs. Substituting (15) into (14) and expressing the sums as Bessel functions, we arrive at the central result of this manuscript, the approximate analytical expression for the glide velocity: The term associated with the formation of kink pairs is not included as it is negligible in the kink-dominated regime z 1. The assumptions leading to the derivation of the drift velocity are only strictly valid in the limit z 1, pertinent to the overwhelming majority of experimental conditions [25]. Leaving out the kink pair formation term has a simple but clear benefit of extending the single kink regime in the drift velocity without the onset of saturation effects at z ∼ 1 and beyond which allows the velocity law to be manually extended at any point by interpolating it to other expressions describing dislocation mobility in the viscous drag or high-speed regime [26]. In the Appendix B we present an expression for glide velocity derived using the exact combinatorics of non-overlapping kinks, and demonstrate that the approximate expression (16) is entirely consistent with the exact approach in the kink-dominated mobility regime. Comparison to kMC.-We implemented a rejectionfree kMC algorithm [31] to propagate the line and extract the mean glide velocity. While kMC does not provide closed form analytical expressions for the expectation values, it serves as a reference for the approximate analytical solution (16). We used the following model rates for the fundamental reactions: and an expression for the free energy of a kink from Ref. [13] A 1 2 111 screw dislocation in iron is described by the following set of atomistically obtained parameters b = 2.47Å, U k = 0.33 eV [32], T ath = 700 K [33], σ P = 900 MPa [34], and kink dissipation γ = 1.83 m u ps −1 [10]. Rates k ± m are derived by equating the mean velocity of a kink under shear stress [10]  k ± m ∼ k 0 (1 ± ασ), and subsequently solving for α. Estimating k 0 from the Debye frequency ∼ 10 ps −1 gives k ± m = k 0 (1 ± 0.004 MPa −1 σ). We assume that the energy barrier to kink pair formation is equal to the free energy of formation, which is consistent with the minimum energy pathways for kink formation obtained from atomistic simulations [22,32,35,36].
We ran kMC simulations for two dislocation segment lengths n = 50 and 250 for a range of stresses and temperatures for which f k > 0, see Fig. 3. Statistical uncertainties in the velocities obtained by kMC simulations are below 3 %. The analytically computed velocities (16) are in quantitative agreement with velocities derived from kMC simulations. An estimate for the critical temperature for a gliding screw dislocation is obtained by analogy with the equilibrium case (8) using the condition n = ξ/ √z * . The largest discrepancy is found for σ = 150 MPa at low temperature, where the analytical expression for velocity overestimates the velocity found in kMC simulations by about 50 %, which is expected since the linear response approximation underlying the derivation assumes the low stress limit.
Using the analytical expression (16), it is possible to predict the velocity of a gliding screw dislocation for any bcc material as a function of temperature, stress, and length; provided a model for the free energy of formation of a kink is available. Parameters U k and σ P are read-   (19). Athermal temperatures and kink mobilities are estimated by scaling the value for iron by the ratio of melting temperatures (T ath = 700 K × T mat melt /T Fe melt ) and shear moduli (α = 0.004 MPa −1 × µ mat /µ Fe ), respectively. ily available from experimental studies [11,25], while we infer T ath and k ± m by scaling the parameters for iron, see Tab. I. The glide velocities of a screw dislocation segment under a resolved shear stress of 50 MPa in bcc metals W, Mo, Nb, V, and Fe, are shown in Fig. 4, assuming attempt rate k 0 = 10 ps −1 and kink height h = 2/3a.
Concluding remarks.-We derived an analytical expression for the thermally activated glide velocity of a screw dislocation segment of arbitrary length in a bcc metal, driven by external stress. The model accounts for formation energies and rates of fundamental reactions involving kinks, which are accessible to atomistic simulations [10,[40][41][42], and as such is generally transferable across bcc metals. We show that the activation energy for dislocation glide is halved at high tempera-ture, as the system crosses a critical transition into the thermodynamic limit. The mobility law generalizes past formulations of kink-limited motion of screw dislocations. We note that dislocation mobility enters the viscous drag regime [26,43] for f k 0, which lies outside the scope of this work.
Our model shows that the activation energy for plastic deformation can indeed be lowered to f k , however, it also suggests that the activation energy should be closer to 2f k in irradiated microstructures, as the mean dislocation length is reduced by radiation defects acting as dislocation pinning points [13]. The dynamics of screw dislocations terminating at features often found in realistic microstructure, such as junctions, grain boundaries, or surfaces, is a largely unexplored area. We hence emphasize the need for atomistic studies of kinetics of such systems in order to fully classify the rate-limiting mechanisms contributing to plastic deformation of bcc metals.  Since in the kink-dominated regime we have √ z 1, the logarithmic term in the above equation can be simplified using a Taylor expansion and neglecting the terms quadratic in √ z, leading to the expression where I 0 (x) is the modified Bessel function of the first kind of order zero. We note that the exact solution in terms of the hypergeometric function is identified as the integral (A4).

Appendix B: Statistics of the exact partition function and comparison to kMC
In the manuscript we defined the critical value T * as the point where the change of slope of ln r is maximal under the approximation that the free energy of formation of a kink pair is given by a linear function of temperature. The same procedure is repeated for the exact partition function (2) with the mean number of kink pairs This time the stationary point depends on system size n, therefore we investigate the equation numerically on the interval n ∈ [10 1 , 10 5 ], leading to the following implicit equation: where ξ(n) in the limit n 1 is well approximated by the series ξ(n) = 0.9548 + 1.2938 n − 7.1232 n 2 + . . .
It is remarkable that for large systems, where n → ∞, the scaling function ξ(n) becomes identical to the numerical constant ξ = 0.9548 determined for the approximate partition function. In the exact partition function we permit up to one kink per site, and hence do not include microstates with overlapping kinks. As the system size becomes large, for a given number of kink pairs r, the number of microstates with overlapping kinks becomes negligible compared to the number of microstates with non-overlapping kinks, which reaffirms that the exact and approximate partition functions are equivalent at low kink densities.
The drift velocity for the system described by the exact partition function is determined in the linear response approximation v eq = Z −1 where Z = 2 F 1 1 2 − n 2 , − n 2 ; 1; 4z is the exact partition function. As before, we put together all the microstates with the same number of kink pairs into a group P r .
It remains to determine the group-averaged line velocity v r (13), which requires knowledge of the total number of microstates connected to each group n ± r,ε = i∈Pr n ± i,ε for each fundamental reaction or process. This combinatorial problem is exactly solvable for non-overlapping kinks. We begin by considering the process of kink pair annihilation, here explicitly included for completeness with a rate consistent with kink motion k ± a = k ± m . Any microstate C i with non-overlapping kinks can equivalently be described by a sequence of horizontal steps h, left kink steps l, and right kink steps r, with the number of l and r steps each being equal to r. A process of kink pair annihilation advancing the line towards the biasing direction is described by a subsequence (l,r) in position (i, i + 1) modulo n turning into the subsequence (h,h), with the rest of the line remaining unchanged. There are n ways to choose i, and of the remaining (n − 2) sites we need to choose (r − 1) steps of l and r each, giving the total number of annihilation processes in P r as where we used n + i,ε = n − i,ε due to symmetry. Next, we recognize that the total number of pair annihilation processes in group P r corresponds precisely to the total number of pair creation processes from P r−1 , leading to n ± r,a = n ∓ r−1,c . We use a similar argument to determine the number of processes of kink motion advancing the line towards the biasing direction, with a process described by either subsequence (r, h) turning into (h, r), or (h, l) turning into (l, h). The resulting group-averaged veloc- where N r = n r n−r r is the number of microstates containing r kink pairs. Substituting (B6) into (B4), we arrive at the analytical expression for the glide velocity: v eq = 2 F 1 1 2 − n 2 , − n 2 ; 1; 4z The above expression is exact in the sense that the combinatorics of non-overlapping kinks is explicitly accounted for.
The derivation of the drift velocity (B7) relies on the linear response approximation and the assumption that all microstates with the same number of kink pairs appear with equal probability. In contrast to the solution presented in the manuscript which is valid for z 1, here we did not assume a low kink density with the caveat that kinks were not permitted to overlap. Consequently, at high temperature the line saturates with kinks, leading to an eventual suppression of kink creation and motion processes. In Fig. 5 we compare the analytical expression for drift velocity (B7) to numerical simulation. Note that the kinetic Monte Carlo simulation here is restricted to only permit one kink per site, in accordance with the combinatorics of non-overlapping kinks. The onset of a saturation in velocity is apparent at high temperature, while at low temperature, where z 1, the velocities are identical to the approximate solution (16).